#include "assert.h"
#include "ctype.h"
#include "float.h"
#include "math.h"
#include "stdio.h"
#include "string.h"
#include "stdlib.h"
#include "time.h"
#include "algorithm"
#include "numeric"
#include "functional"
#include "utility"
#include "bitset"
#include "vector"
#include "list"
#include "set"
#include "map"
#include "queue"
#include "stack"
#include "string"
#include "sstream"
#include "iostream"
using namespace std;

#pragma warning(disable:4018)  // signed/unsigned mistatch
#pragma warning(disable:4244)  // w64 to int cast
#pragma warning(disable:4267)  // big to small -- possible loss of data
#pragma warning(disable:4786)  // long identifiers
#pragma warning(disable:4800)  // forcing int to bool
#pragma warning(disable:4996)  // deprecations

typedef long long i64;
#define all(x) (x).begin(), (x).end()

namespace acm {

///////////////////// parsing functions ////////////////////////////

char* gstrsep(char** nptr, const char* delim) {
  while (**nptr != '\0' && strchr(delim, **nptr) != NULL) ++(*nptr);
  if (**nptr == '\0') return NULL;
  char* const token = (*nptr)++;
  while (**nptr != '\0' && strchr(delim, **nptr) == NULL) ++(*nptr);
  if (**nptr != '\0') *(*nptr)++ = '\0';
  return token;
}

////////////////////// floating point /////////////////////////////

const double eps = 5e-8;
int DoubleComp(double a, double b) {
  if (fabs(a - b) < eps) return 0;
  if (b * (1.0 - eps) < a && a < b * (1.0 + eps)) return 0;
  if (b * (1.0 + eps) < a && a < b * (1.0 - eps)) return 0;
  return a < b ? -1 : +1;
}

int sign(int x) {return x < 0 ? -1 : +1;}
int DoubleSign(double x) {return DoubleComp(x, 0.0) < 0 ? -1 : +1;}
int DoubleCast(double x) {return (int)(x + (x < 0.0 ? -eps : +eps));}
int DoubleFloor(double x) {return DoubleCast(floor(x + eps));}
int DoubleCeil(double x) {return DoubleCast(ceil(x - eps));}

///////////////////// date-related ////////////////////////////

static const int kMonthDays[12] = {
  31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31
};
static const int kMonthNumber[2][12] = {
  {0, 3, 3, 6, 1, 4, 6, 2, 5, 0, 3, 5},
  {6, 2, 3, 6, 1, 4, 6, 2, 5, 0, 3, 5}
};

bool IsLeapYear(int year) {
  return year % 4 == 0 && (year % 100 != 0 || year % 400 == 0);
}

int DayOfTheWeek(int year, int month, int day) {
  // returns the week day of the specified date, where 0 is Monday,
  // 1 is Tuesday and so on...
  const int isLeap = IsLeapYear(year) ? 1 : 0;
  const int ncentury = (3 - (year / 100) % 4) * 2;
  const int nyear = 5 * (year % 100) / 4;
  const int nmonth = kMonthNumber[isLeap][month - 1];
  return (ncentury + nyear + nmonth + day - 1) % 7;
}

///////////////////// infix to postfix ////////////////////////////

string InfixToPostfix(const string& infix) {
  stack<char> operand; string postfix;
  for (int k = 0; k < infix.size(); ++k) {
    if (infix[k] == '(') {
      operand.push('(');
    } else if (infix[k] == ')') {
      while (!operand.empty() && operand.top() != '(') {
        postfix.push_back(operand.top()); operand.pop();
      }
      operand.pop();
    } else if (strchr("*/", infix[k]) != NULL) {
      while (!operand.empty() && strchr("+-(", operand.top()) == NULL) {
        postfix.push_back(operand.top()); operand.pop();
      }
      operand.push(infix[k]);
    } else if (strchr("+-", infix[k]) != NULL) {
      while (!operand.empty() && operand.top() != '(') {
        postfix.push_back(operand.top()); operand.pop();
      }
      operand.push(infix[k]);
    } else {
      postfix.push_back(infix[k]);
    }
  }
  while (!operand.empty()) {
    if (operand.top() != '(') postfix.push_back(operand.top());
    operand.pop();
  }
  return postfix;
}

///////////////////// knuth - morris - pratt ////////////////////////////

void ConstructFailTable(const string& W, vector<int>& failTable) {
  const int N = W.size(); failTable.assign(N, 0);
  failTable[0] = -1, failTable[1] = 0;
  for (int k = 2, i = 0; k < N;) {
    if (W[k - 1] == W[i]) failTable[k++] = ++i;
    else if (i > 0) i = failTable[i];
    else failTable[i++] = 0;
  }
}

int SearchWord(const string& S, const string& W, const vector<int>& failTable) {
  for (int index = 0, k = 0; index + k < S.size();) {
    if (k >= W.size()) return index;
    if (W[k] == S[index + k]) ++k;
    else {
      index += k - failTable[k];
      if (k > 0) k = failTable[k];
    }
  }
  return string::npos;
}

///////////////////// aho - corasick ////////////////////////////

static const int MAX_PATTERN_LENGTH = 1001;  // maximum pattern length
static const int SYMBOL_COUNT = 52;          // number of possible symbols

struct node_t {
  static const int ROOT = 0;
  node_t() : failure(-1) {memset(next, -1, sizeof(next));}
  int failure;             // failure node
  int next[SYMBOL_COUNT];  // node out-links
  set<int> output;         // list of patterns found
};
typedef vector<node_t> automaton_t;

int GetSymbolValue(char c) {
  if (isupper(c)) return c - 'A';
  else return c - 'a' + ('Z' - 'A' + 1);
}

void CreateKeywordTree(automaton_t& automaton, const char pattern[][MAX_PATTERN_LENGTH], int N) {
  automaton.assign(1, node_t());
  for (int k = 0; k < N; ++k) {
    int node = node_t::ROOT;
    for (const char* text = pattern[k]; *text != '\0'; ++text) {
      const int index = GetSymbolValue(*text);
      if (automaton[node].next[index] != -1) {
        node = automaton[node].next[index];
      } else {
        const int newNode = (int)automaton.size();
        automaton.push_back(node_t());
        automaton[node].next[index] = node = newNode;
      }
    }
    automaton[node].output.insert(k);
  }
}

void CreateFailureLinks(automaton_t& automaton) {
  queue<int> nqueue; nqueue.push(node_t::ROOT);
  while (!nqueue.empty()) {
    const int node = nqueue.front(); nqueue.pop();
    for (int k = 0; k < SYMBOL_COUNT; ++k) {
      const int next = automaton[node].next[k];
      if (next != -1) {
        nqueue.push(next);
        int failure = automaton[node].failure;
        while (failure != -1 && automaton[failure].next[k] == -1)
          failure = automaton[failure].failure;
        automaton[next].failure = failure == -1 ? node_t::ROOT : automaton[failure].next[k];
        automaton[next].output.insert(all(automaton[automaton[next].failure].output));
      }
    }
  }
}

void CreateAutomaton(automaton_t& automaton, const char pattern[][MAX_PATTERN_LENGTH], int N) {
  CreateKeywordTree(automaton, pattern, N);
  CreateFailureLinks(automaton);
}

void SearchPatterns(automaton_t& automaton, const char* str, set<int>& found) {
  int state = node_t::ROOT;
  for (; *str != '\0'; ++str) {
    const int index = GetSymbolValue(*str);
    while (automaton[state].next[index] == -1 && automaton[state].failure != -1)
      state = automaton[state].failure;
    if (automaton[state].next[index] != -1)
      state = automaton[state].next[index];
    found.insert(all(automaton[state].output));
  }
}

///////////////////// binary & ternary search ////////////////////////////

template <class T, class Iterator, class Pred>
Iterator BinarySearch(const T& key, Iterator begin, Iterator end, Pred pred) {
  while (begin < end) {
    const Iterator middle = begin + (end - begin) / 2;
    if (pred(*middle, key)) begin = middle + 1;
    else end = middle;
  }
  return begin;
}

template <class T, class Function>
T TernarySearch(T nmin, T nmax, const Function& func, double eps = 5e-8) {
  for (int k = 0; k < 300 && nmax - nmin >= eps; ++k) {
    const T nleft = (nmin * 2 + nmax) / 3;
    const T nright = (nmin + nmax * 2) / 3;
    if (func(nleft) < func(nright)) nmin = nleft;
    else nmax = nright;
  }
  return (nmax + nmin) / 2;
}

///////////////////// bubble sort ////////////////////////////

template <class Iterator, class Pred>
int BubbleSort(Iterator begin, Iterator end, Pred pred) {
  int swapCount = 0; bool changed = true;
  for (Iterator nend = --end; changed && begin != nend; --nend) {
    changed = false;
    for (Iterator nptr = begin; nptr != nend; ++nptr) {
      if (!pred(*nptr, *(nptr + 1)))
        swap(*nptr, *(nptr + 1)), ++swapCount, changed = true;
    }
  }
  return swapCount;
}

////////////////////// merge sort ////////////////////////////

template <class T, class Iterator, class Pred>
int Merge(Iterator begin, Iterator middle, Iterator end, const Pred& pred, vector<T>& aux) {
  int swapCount = 0; aux.clear();
  Iterator left = begin, right = middle;
  while (left != middle && right != end) {
    if (pred(*left, *right)) aux.push_back(*left++);
    else aux.push_back(*right++), swapCount += middle - left;
  }
  copy(left, middle, back_inserter(aux));
  copy(right, end, back_inserter(aux));
  copy(aux.begin(), aux.end(), begin);
  return swapCount;
}

template <class T, class Iterator, class Pred>
int MergeSort(Iterator begin, Iterator end, const Pred& pred, vector<T>& aux) {
  const int size = end - begin;
  if (size <= 1) return 0;
  const Iterator middle = begin + size / 2;
  int swapCount = MergeSort<T, Iterator, Pred>(begin, middle, pred, aux);
  swapCount += MergeSort<T, Iterator, Pred>(middle, end, pred, aux);
  swapCount += Merge<T, Iterator, Pred>(begin, middle, end, pred, aux);
  return swapCount;
}

template <class T, class Iterator, class Pred>
int MergeSort(Iterator begin, Iterator end, Pred pred) {
  vector<T> aux; aux.reserve(end - begin);
  return MergeSort(begin, end, pred, aux);
}

////////////////////// Range Minimum Query (RMQ) ////////////////////////////

namespace rmq {

const int MAX_N = 100000, MAX_LOG_N = 20;
template <class T> void PreprocessRMQ(const vector<T>& V, int M[MAX_N][MAX_LOG_N]) {
  // creates a table where the position (i,j) have the index of the minimum
  // element in V in the interval [i,i+2^j)
  const int N = V.size();
  for (int i = 0, N = V.size(); i < N; ++i)
    M[i][0] = i;
  for (int j = 1; (1 << j) <= N; ++j) for (int i = 0; i + (1 << j) <= N; ++i) {
    const int left = M[i][j - 1], right = M[i + (1 << (j - 1))][j - 1];
    M[i][j] = V[left] <= V[right] ? left : right;
  }
}

template <class T> int QueryRMQ(const vector<T>& V, const int M[MAX_N][MAX_LOG_N], int i, int j) {
  // retrieves the index of the minimum element in V in the interval [i,j)
  int k = 0; for (int d = j - i; d > 1; d >>= 1, ++k);
  const int left = M[i][k], right = M[j - (1 << k)][k];
  return V[left] <= V[right] ? left : right;
}

};

////////////////////// big mod & power ////////////////////////////

template <class T> T BigMod(T b, T p, T m) {
  if (m == 1) return 0;
  if (p == 0) return 1;
  if (p == 1) return b % m;
  T x = BigMod(b, p / 2, m);
  x = (x * x) % m;
  if (p % 2 != 0) x = (x * (b % m)) % m;
  return x;
}

template <class T> T Power(T b, T p) {
  if (p == 0) return 1;
  if (p == 1) return b;
  T x = Power(b, p / 2); x = x * x;
  if (p % 2 != 0) x *= b; return x;
}

////////////////////// GCD ////////////////////////////

template <class T> T GCD(T a, T b) {
  if (a < 0) a = -a; if (b < 0) b = -b;
  if (b > a) swap(a, b);
  while (b != 0) {a %= b; swap(a, b);}
  return a == 0 ? 1 : a;
}

template <class T> T LCM(T a, T b) {
  return (a / GCD(a, b)) * b;
}

template <class T> T DivByGCD(T& a, T& b) {
  const T d = GCD(a, b);
  return a /= d, b /= d, d;
}

template <class T> T Modulo(T x, T n) {
  return (n + x % n) % n;
}

////////////////////// modular inverse ////////////////////////////

template <class T> T ExtendedGCD(T a, T b, T& x, T& y) {
  // xa + yb = GCD(a, b)
  if (b > a) return ExtendedGCD(b, a, y, x);
  if (b == 0) return x = 1, y = 0, a;
  T x1, y1, d = ExtendedGCD(b, a % b, x1, y1);
  return x = y1, y = x1 - (a / b) * y1, d;
}

template <class T> pair<T, bool> ModularInverse(T a, T b, T n) {
  // ax = b (mod n), all solutions are of the form x + k * (n / GCD(a, n))
  a = Modulo(a, n), b = Modulo(b, n);
  T x, y, d = ExtendedGCD(a, n, x, y);
  if (b % d != 0) return make_pair(0, false);
  else return make_pair(Modulo((x * b) / d, n), true);
}

////////////////////// chinese remainder theorem ////////////////////////////

template <class T> bool ChineseRemainder(T a1, T n1, T a2, T n2, T& a3, T& n3) {
  T r, s, d = ExtendedGCD(n1, n2, r, s);
  a3 = Modulo(a1 * s * (n2 / d) + a2 * r * (n1 / d), n3 = (n1 / d) * n2);
  return a1 % d == a2 % d;
}

template <class T> bool ChineseRemainder(const vector<T>& A, const vector<T>& N, T& a0, T& n0) {
  T a1, n1; a0 = 0, n0 = 1;
  for (int k = 0; k < N.size(); ++k) {
    if (ChineseRemainder(A[k], N[k], a0, n0, a1, n1)) a0 = a1, n0 = n1;
    else return false;
  }
  return true;
}

////////////////////// combinatorial recurrences ////////////////////////////

template <class T> T IterativeChoose(T n, T k) {
  // how many sets of K objects can be taken from a group of N different
  // objects, where order doesn't matter and objects aren't replaced
  if (k > n) return 0; if (k > n / 2) k = n - k;
  T num = 1, den = 1;
  for (T i = k; i > 0; --i) {
    T toMul = n - k + i, toDiv = i;
    DivByGCD(toMul, toDiv); DivByGCD(num, toDiv); DivByGCD(toMul, den);
    num *= toMul, den *= toDiv;
  }
  return num / den;
}

static const int MAX_N = 50;
static long long choose[MAX_N][MAX_N];  // <-- initialize with -1
long long Choose(int N, int K) {
  // how many sets of K objects can be taken from a group of N different
  // objects, where order doesn't matter and objects aren't replaced
  if (K == 0 || K == N) return 1;
  if (choose[N][K] != -1) return choose[N][K];
  return choose[N][K] = Choose(N - 1, K) + Choose(N - 1, K - 1);
}

long long MultiChoose(int N, int K) {
  // how many sets of K objects can be taken from a group of N different
  // objects, where order doesn't matter and objects are replaced
  return Choose(N + K - 1, K);
}

long long OrderedChoose(int N, int K) {
  // how many sets of K objects can be taken from a group of N different
  // objects, where order matters and objects aren't replaced
  long long count = 1;
  for (int k = 0; k < K; ++k) count *= (N - k);
  return count;
}

long long MultiOrderedChoose(int N, int K) {
  // how many sets of K objects can be taken from a group of N different
  // objects, where order matters and objects are replaced
  return Power(N, K);
}

static long long stirling[MAX_N][MAX_N];  // <-- initialize with -1
long long Stirling(int N, int K) {
  // how many different K groups can be created from N different elements (order doesn't matter)
  if (K == 1 || K == N) return 1;
  if (stirling[N][K] != -1) return stirling[N][K];
  return stirling[N][K] = K * Stirling(N - 1, K) + Stirling(N - 1, K - 1);
}

static long long euler[MAX_N][MAX_N];  // <-- initialize with -1
long long Euler(int N, int K) {
  // how many permutations of the sequence {0, 1, ..., N - 1} have
  // K ascents (i.e.: positions where S[k] < S[k + 1])
  if (K == 0 || K == N - 1) return 1;
  if (euler[N][K] != -1) return euler[N][K];
  return euler[N][K] = (N - K) * Euler(N - 1, K - 1) + (K + 1) * Euler(N - 1, K);
}

static long long derangment[MAX_N];  // <-- initialize with -1
long long Derangment(int N) {
  // how many permutations of N elements exists such that for every k in 1..N
  // P[k] != k, i.e.: there are no fixed points or cycles of length 1
  if (N == 0 || N == 1) return 1 - N;
  if (derangment[N] != -1) return derangment[N];
  return derangment[N] = (N - 1) * (Derangment(N - 1) + Derangment(N - 2));
}

static long long mergeCount[MAX_N][MAX_N];  // <-- initialize with -1
long long MergeCount(int N, int K) {
  if (N == 0 || K == 0) return 1;
  if (mergeCount[N][K] != -1) return mergeCount[N][K];
  mergeCount[N][K] = 0;
  for (int i = 0; i <= N; ++i)
    mergeCount[N][K] += MergeCount(N - i, K - 1);
  return mergeCount[N][K];
}

long long Catalan(int N) {
  static const int MAX_N = 10000;  // max value to calculate
  static long long catalan[MAX_N] = {1, 1}; static int i = 2;
  for (; i <= N; ++i) catalan[i] = (catalan[i - 1] * (4 * i - 2)) / (i + 1);
  return catalan[N];
}

////////////////////// prime numbers ////////////////////////////

void GeneratePrimes(int N, vector<bool>& primeMark, vector<int>& primes) {
  primeMark.clear(), primeMark.resize(N, true), primeMark[1] = false;
  for (int i = 4; i < N; i += 2) primeMark[i] = false; primes.push_back(2);
  for (long long i = 3; i < N; i += 2) if (primeMark[i]) {
    for (long long j = i * i; j < N; j += 2 * i) primeMark[j] = false;
    primes.push_back((int)i);
  }
}

typedef pair<int, int> factor_t;  // (factor, count)
void PrimeFactors(int N, const vector<int>& primes, vector<factor_t>& factors) {
  const int sq = DoubleFloor(sqrt((double)N));
  for (int k = 0; k < primes.size() && primes[k] <= sq && N > 1; ++k) {
    int count = 0;
    for (; N % primes[k] == 0; ++count, N /= primes[k]);
    if (count > 0) factors.push_back(factor_t(primes[k], count));
  }
  if (N > 1) factors.push_back(factor_t(N, 1));
}

int GetFactorCount(int N, int P) {
  int count = 0;
  for (int h = P; h <= N; h *= P) count += N / h;
  return count;
}

void FactorialFactors(int N, const vector<int>& primes, vector<factor_t>& factors) {
  for (int k = 0; k < primes.size() && primes[k] <= N; ++k)
    factors.push_back(factor_t(primes[k], GetFactorCount(N, primes[k])));
}

int EulerPhi(int N, const vector<int>& primes, vector<factor_t>& factors) {
  // how many numbers less than N are coprime with N
  PrimeFactors(N, primes, factors);
  int coprimeCount = 1;
  for (int k = 0; k < factors.size(); ++k) {
    const int prime = factors[k].first, exp = factors[k].second;
    coprimeCount *= Power(prime, exp - 1) * (prime - 1);
  }
  return coprimeCount;
}

////////////////////// matrix and vector operations ////////////////////////////

namespace calculus {

typedef vector<double> vector_t;
typedef vector<vector_t> matrix_t;

matrix_t IdentityMatrix(int N) {
  matrix_t I(N, vector_t(N, 0));
  for (int k = 0; k < N; ++k) I[k][k] = 1;
  return I;
}

vector_t MultVector(const matrix_t& A, const vector_t& X) {
  // A * X = B
  const int M = A.size(), N = X.size();
  vector_t B(M, 0);
  for (int i = 0; i < M; ++i) for (int k = 0; k < N; ++k)
    B[i] += A[i][k] * X[k];
  return B;
}

matrix_t MultMatrix(const matrix_t& A, const matrix_t& B) {
  // A * B = C
  const int M = A.size(), L = B.size(), N = B[0].size();
  matrix_t C(M, vector_t(N, 0));
  for (int i = 0; i < M; ++i) for (int j = 0; j < N; ++j)
    for (int k = 0; k < L; ++k) C[i][j] += A[i][k] * B[k][j];
  return C;
}

matrix_t PowerMatrix(const matrix_t& B, int P) {
  // B ^ P = A
  const int N = B.size();
  if (P == 0) return IdentityMatrix(N);
  if (P == 1) return B;
  matrix_t A = PowerMatrix(B, P / 2);
  A = MultMatrix(A, A);
  if (P % 2 == 1)
    A = MultMatrix(A, B);
  return A;
}

////////////////////// linear equation systems ////////////////////////////

void FactorizationLU(const matrix_t& A, matrix_t& L, matrix_t& U) {
  const int N = A.size(); L.assign(N, vector_t(N, 0.0)); U = A;
  for (int k = 0; k < N; ++k) L[k][k] = 1.0;
  for (int i1 = 0; i1 < N; ++i1) for (int i2 = i1 + 1; i2 < N; ++i2) {
    L[i2][i1] = U[i2][i1] / U[i1][i1];
    for (int j = i1; j < N; ++j)
      U[i2][j] -= L[i2][i1] * U[i1][j];
  }
}

void SolveSystem(const matrix_t& L, const matrix_t& U, const vector_t& B, vector_t& X) {
  const int N = L.size(); vector_t Y(N, 0.0);
  for (int i = 0; i < N; ++i) {
    Y[i] = B[i];
    for (int j = 0; j < i; ++j)
      Y[i] -= L[i][j] * Y[j];
    Y[i] /= L[i][i];
  }
  X.assign(N, 0.0);
  for (int i = N - 1; i >= 0; --i) {
    X[i] = Y[i];
    for (int j = i + 1; j < N; ++j)
      X[i] -= U[i][j] * X[j];
    X[i] /= U[i][i];
  }
}

////////////////////// matrix determinant ////////////////////////////

double GetDeterminant(matrix_t& M) {
  const int N = M.size(); double determinant = 1.0;
  for (int i = 0; i < N; ++i) {
    if (M[i][i] == 0) {
      int row = -1;
      for (int ni = i + 1; ni < N && row == -1; ++ni)
        if (M[ni][i] != 0) row = ni;
      if (row == -1) return 0.0;
      swap(M[i], M[row]);
      determinant *= -1;
    }
    for (int ni = i + 1; ni < N; ++ni) {
      const double Lij = M[ni][i] / M[i][i];
      for (int j = i; j < N; ++j)
        M[ni][j] -= M[i][j] * Lij;
    }
  }
  for (int i = 0; i < N; ++i)
    determinant *= M[i][i];
  return determinant;
}

////////////////////// two-phase simplex ////////////////////////////

template <class T> struct var_t {
  var_t(T value_ = T(0), bool basic_ = false)
    : value(value_), basic(basic_) {}
  T value; bool basic;
};
template <class T> struct table_t {
  static const int MAX_M = 100;  // max constraints
  static const int MAX_N = 100;  // max variables
  T t[MAX_M][MAX_N]; int M, N;
};

template <class T>
int GetInColumn(const table_t<T>& table, int M, int N) {
  int col = -1;
  for (int j = 0; j < N; ++j) if (table.t[M][j] < T(0))
    if (col == -1 || table.t[M][j] < table.t[M][col]) col = j;
  return col;
}

template <class T>
int GetOutRow(const table_t<T>& table, int M, int N, int col) {
  int row = -1; T minCoef = T(0);
  for (int i = 0; i < M; ++i) if (table.t[i][col] > T(0)) {
    const T coef = table.t[i][N] / table.t[i][col];
    if (row == -1 || coef < minCoef)
      row = i, minCoef = coef;
  }
  return row;
}

template <class T>
void MakeGauss(table_t<T>& table, int M, int N, int row, int col) {
  const T pivot = table.t[row][col];
  for (int j = 0; j <= N; ++j) table.t[row][j] /= pivot;
  for (int i = 0; i <= M; ++i) if (i != row) {
    const T factor = table.t[i][col];
    for (int j = 0; j <= N; ++j)
      table.t[i][j] -= table.t[row][j] * factor;
  }
}

template <class T>
var_t<T> GetVariable(const table_t<T>& table, int M, int N, int col) {
  bool ok = true; int row = -1;
  for (int i = 0; i <= M && ok; ++i) if (table.t[i][col] != T(0))
    if (row == -1) row = i; else ok = false;
  if (row == -1 || !ok) return var_t<T>(T(0), false);
  else return var_t<T>(table.t[row][N] / table.t[row][col], true);
}

template <class T>
bool Simplex(table_t<T>& table, int M, int N, int C, var_t<T>* vars) {
  for (int col; (col = GetInColumn(table, M, N)) != -1;) {
    const int row = GetOutRow(table, C, N, col);
    if (row == -1) return false;
    MakeGauss(table, M, N, row, col);
  }
  for (int j = 0; j < N; ++j)
    vars[j] = GetVariable(table, M, N, j);
  vars[N] = var_t<T>(table.t[M][N], true);
  return true;
}

template <class T>
bool Simplex(table_t<T>& table, var_t<T>* vars) {
  return Simplex(table, table.M, table.N, table.M, vars);
}

template <class T>
bool TwoPhaseSimplex(table_t<T>& table, var_t<T>* vars) {
  for (int j = 0; j <= table.N; ++j) {
    table.t[table.M + 1][j] = T(0);
    for (int i = 0; i < table.M; ++i)
      table.t[table.M + 1][j] -= table.t[i][j];
  }
  bool ok = Simplex(table, table.M + 1, table.N, table.M, vars);
  if (ok) ok = table.t[table.M + 1][table.N] == T(0);
  return ok && Simplex(table, vars);
}

};

////////////////////// disjoint sets ////////////////////////////

int GetNodeRoot(vector<int>& parent, int node) {
  if (parent[node] == -1) return node;
  return parent[node] = GetNodeRoot(parent, parent[node]);
}

bool JoinNodes(vector<int>& parent, int node1, int node2) {
  const int root1 = GetNodeRoot(parent, node1);
  const int root2 = GetNodeRoot(parent, node2);
  if (root1 == root2) return false;
  return parent[root1] = root2, true;
}

////////////////////// depth first search ////////////////////////////

namespace dfs {

struct edge_t {
  edge_t(int src_, int dest_)
    : src(src_), dest(dest_) {}
  int src, dest;
};
typedef vector<edge_t> node_t;
struct graph_t {
  static const int MAX_N = 500;
  node_t V[MAX_N]; int N;
};

struct dfs_t {
  dfs_t(int N) : parent(N, -1), preindex(N, -1), postindex(N, -1), lowindex(N, N) {}
  vector<int> parent, preorder, preindex, postorder, postindex, lowindex;
};

void DepthFirstSearch(const graph_t& graph, dfs_t& dfs, int parent, int node) {
  dfs.parent[node] = parent;
  dfs.preindex[node] = (int)dfs.preorder.size();
  dfs.preorder.push_back(node);
  for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
    const edge_t* const edge = &graph.V[node][k];
    if (dfs.parent[edge->dest] == -1) { /* on forward edge... */
      DepthFirstSearch(graph, dfs, node, edge->dest);
    } else {
      if (dfs.postindex[edge->dest] == -1) /* on back edge... */;
      else if (dfs.preindex[edge->dest] > dfs.preindex[node]) /* on down edge... */;
      else /* on cross edge... */;
    }
    dfs.lowindex[node] = min(dfs.lowindex[node], dfs.lowindex[edge->dest]);
  }
  dfs.postindex[node] = (int)dfs.postorder.size();
  dfs.postorder.push_back(node);
}

void DepthFirstSearch(const graph_t& graph, dfs_t& dfs) {
  for (int k = 0; k < graph.N; ++k) if (dfs.parent[k] == -1)
    DepthFirstSearch(graph, dfs, k, k);
}

////////////////////// kosaraju's strongly connected components ////////////////////////////

struct kernel_t {
  static const int MAX_N = graph_t::MAX_N;
  kernel_t(int N) : componentIndex(N, -1) {}
  graph_t graph;
  vector<int> componentIndex;
  vector<int> componentNodes[MAX_N];
};

void KosarajuDFS(const graph_t& graph, kernel_t& kernel, dfs_t& dfs, int component, int parent, int node) {
  kernel.componentNodes[component].push_back(node);
  kernel.componentIndex[node] = component;
  dfs.parent[node] = parent;
  dfs.preindex[node] = (int)dfs.preorder.size();
  dfs.preorder.push_back(node);
  for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
    const edge_t* const edge = &graph.V[node][k];
    if (dfs.parent[edge->dest] == -1) {
      KosarajuDFS(graph, kernel, dfs, component, node, edge->dest);
    } else if (dfs.postindex[edge->dest] != -1 && dfs.preindex[edge->dest] <= dfs.preindex[node]) {
      const int destComponent = kernel.componentIndex[edge->dest];
      if (component != destComponent)
        kernel.graph.V[component].push_back(edge_t(component, destComponent));
    }
  }
  dfs.postindex[node] = (int)dfs.postorder.size();
  dfs.postorder.push_back(node);
}

void KosarajuDFS(const graph_t& graph, kernel_t& kernel, dfs_t& rdfs) {
  dfs_t dfs(graph.N);
  for (int k = graph.N - 1; k >= 0; --k) {
    const int node = rdfs.postorder[k];
    if (dfs.parent[node] != -1) continue;
    const int component = kernel.graph.N++;
    KosarajuDFS(graph, kernel, dfs, component, node, node);
  }
}

void KosarajuSCC(const graph_t& graph, const graph_t& rgraph, kernel_t& kernel) {
  kernel.graph.N = 0;
  dfs_t rdfs(rgraph.N);
  DepthFirstSearch(rgraph, rdfs);
  KosarajuDFS(graph, kernel, rdfs);
}

////////////////////// tarjan's strongly connected components ////////////////////////////

struct EdgeFunctor : public binary_function<edge_t, edge_t, bool> {
  bool operator ()(const edge_t& e1, const edge_t& e2) {
    if (e1.src != e2.src) return e1.src < e2.src;
    else return e1.dest < e2.dest;
  }
};

void ConstructKernel(const graph_t& graph, kernel_t& kernel) {
  for (int k = 0; k < kernel.graph.N; ++k) kernel.graph.V[k].clear();
  for (int k = 0; k < graph.N; ++k) {
    const int src = kernel.componentIndex[k];
    for (int i = 0, N = (int)graph.V[k].size(); i < N; ++i) {
      const edge_t* const edge = &graph.V[k][i];
      const int dest = kernel.componentIndex[edge->dest]; if (src == dest) continue;
      const node_t::iterator it = lower_bound(all(kernel.graph.V[src]), edge_t(src, dest), EdgeFunctor());
      if (it == kernel.graph.V[src].end() || it->dest != dest)
        kernel.graph.V[src].insert(it, edge_t(src, dest));
    }
  }
}

void TarjanDFS(const graph_t& graph, kernel_t& kernel, dfs_t& dfs, stack<int>& path, int node) {
  dfs.lowindex[node] = dfs.preindex[node] = (int)dfs.preorder.size();
  dfs.preorder.push_back(node); path.push(node);
  for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
    const edge_t* const edge = &graph.V[node][k];
    if (dfs.preindex[edge->dest] == -1)
      TarjanDFS(graph, kernel, dfs, path, edge->dest);
    if (kernel.componentIndex[edge->dest] == -1)
      dfs.lowindex[node] = min(dfs.lowindex[node], dfs.lowindex[edge->dest]);
  }
  if (dfs.lowindex[node] == dfs.preindex[node]) {
    const int component = kernel.graph.N++;
    while (true) {
      const int nextNode = path.top(); path.pop();
      kernel.componentIndex[nextNode] = component;
      kernel.componentNodes[component].push_back(nextNode);
      if (nextNode == node) break;
    }
  }
}

void TarjanSCC(const graph_t& graph, kernel_t& kernel) {
  stack<int> path; kernel.graph.N = 0; dfs_t dfs(graph.N);
  for (int k = 0; k < graph.N; ++k) if (dfs.preindex[k] == -1)
    TarjanDFS(graph, kernel, dfs, path, k);
}

////////////////////// gabow's strongly connected components ////////////////////////////

void GabowDFS(const graph_t& graph, kernel_t& kernel, dfs_t& dfs,
              stack<int>& roots, stack<int>& path, int node) {
  dfs.preindex[node] = (int)dfs.preorder.size();
  dfs.preorder.push_back(node);
  path.push(node); roots.push(node);
  for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
    const edge_t* const edge = &graph.V[node][k];
    if (dfs.preindex[edge->dest] == -1)
      GabowDFS(graph, kernel, dfs, roots, path, edge->dest);
    else if (kernel.componentIndex[edge->dest] == -1)
      while (dfs.preindex[roots.top()] > dfs.preindex[edge->dest]) roots.pop();
  }
  if (roots.top() == node) {
    const int component = kernel.graph.N++; roots.pop();
    while (true) {
      const int nextNode = path.top(); path.pop();
      kernel.componentIndex[nextNode] = component;
      kernel.componentNodes[component].push_back(nextNode);
      if (nextNode == node) break;
    }
  }
}

void GabowSCC(const graph_t& graph, kernel_t& kernel) {
  stack<int> roots, path;
  kernel.graph.N = 0; dfs_t dfs(graph.N);
  for (int k = 0; k < graph.N; ++k) if (dfs.preindex[k] == -1)
    GabowDFS(graph, kernel, dfs, roots, path, k);
}

};

////////////////////// articulations & bridges ////////////////////////////

namespace connectivity {  // <-- pass to static graph

struct DFSData {
  DFSData(int N) : parent(N, -1), preindex(N, -1), lowindex(N, -1) {}
  vector<int> parent, preorder, preindex, lowindex, articulations;
  vector<pair<int, int> > bridges;
};

typedef vector<vector<int> > graph_t;
int DepthFirstSearch(const graph_t& graph, DFSData& dfsData, int parent, int node) {
  dfsData.parent[node] = parent;
  dfsData.lowindex[node] = dfsData.preindex[node] = dfsData.preorder.size();
  dfsData.preorder.push_back(node);
  int biconnectedComponents = 0, subComponentCount = 0;
  bool isArticulation = false;
  for (int k = 0; k < graph[node].size(); ++k) {
    const int dest = graph[node][k];
    if (dfsData.preindex[dest] == -1) {
      subComponentCount++;
      biconnectedComponents += DepthFirstSearch(graph, dfsData, node, dest);
      if (dfsData.lowindex[dest] < dfsData.lowindex[node])
        dfsData.lowindex[node] = dfsData.lowindex[dest];
      if (dfsData.preindex[dest] == dfsData.lowindex[dest])
        dfsData.bridges.push_back(make_pair(node, dest));
      if (dfsData.lowindex[dest] >= dfsData.preindex[node]) biconnectedComponents++;
      if (parent == -1) isArticulation |= subComponentCount >= 2;
      else isArticulation |= dfsData.lowindex[dest] >= dfsData.preindex[node];
    }
    if (dest != parent && dfsData.preindex[dest] < dfsData.lowindex[node])
      dfsData.lowindex[node] = dfsData.preindex[dest];
  }
  if (isArticulation) dfsData.articulations.push_back(node);
  return biconnectedComponents;
}

int DepthFirstSearch(const graph_t& graph, DFSData& dfsData) {
  int biconnectedComponents = 0;
  for (int k = 0; k < graph.size(); k++) {
    if (dfsData.preindex[k] != -1) continue;
    const int subBiconnectedComponents = DepthFirstSearch(graph, dfsData, -1, k);
    biconnectedComponents += max(subBiconnectedComponents, 1);
  }
  return biconnectedComponents;
}

};

////////////////////// topological sorting ////////////////////////////

namespace dfs {

void TopologicalSorting(const graph_t& graph, vector<int>& inDegree, vector<int>& order, vector<int>& index) {
  queue<int> nqueue; index.assign(graph.N, -1); order.clear();
  for (int k = 0, N = (int)inDegree.size(); k < N; ++k)
    if (inDegree[k] == 0) nqueue.push(k);
  for (int k = 0; !nqueue.empty(); ++k) {
    const int node = nqueue.front(); nqueue.pop();
    order.push_back(node), index[node] = k;
    for (int i = 0, N = (int)graph.V[node].size(); i < N; ++i) {
      const edge_t* const edge = &graph.V[node][i];
      if (--inDegree[edge->dest] == 0) nqueue.push(edge->dest);
    }
  }
}

};

////////////////////// kruskal w/ partial sorting ////////////////////////////

namespace kruskal {

struct edge_t {
  edge_t(int u_, int v_, int cost_)
    : u(u_), v(v_), cost(cost_) {}
  int u, v, cost;
};

struct EdgeFunctor : public binary_function<edge_t, edge_t, bool> {
  bool operator ()(const edge_t& e1, const edge_t& e2) const {
    return e1.cost < e2.cost;
  }
};

template <class Iterator>
void JoinAll(Iterator begin, Iterator end, vector<int>& parent, vector<edge_t>& MST) {
  for (; begin != end && MST.size() < parent.size() - 1; ++begin) {
    if (JoinNodes(parent, begin->u, begin->v))
      MST.push_back(*begin);
  }
}

template <class Iterator, class Pred>
void Kruskal(Iterator begin, Iterator end, const Pred& pred, vector<int>& parent, vector<edge_t>& MST) {
  static const int threshold = 100;  // max links to change to normal method
  if (end - begin <= threshold) {
    sort(begin, end, pred);
    JoinAll(begin, end, parent, MST);
  } else {
    const int pivot = begin->cost;
    Iterator middle = partition(begin, end, bind2nd(pred, edge_t(pivot, 0, 0)));
    if (middle == begin) ++middle;
    Kruskal(begin, middle, pred, parent, MST);
    if (MST.size() < parent.size() - 1)
      Kruskal(middle, end, pred, parent, MST);
  }
}

};

////////////////////// prim ////////////////////////////

namespace prim {

const int INFINITE = 100000;
struct edge_t {
  edge_t(int src_, int dest_, int cost_)
    : src(src_), dest(dest_), cost(cost_) {}
  int src, dest, cost;
};
typedef vector<edge_t> node_t;
struct graph_t {
  static const int MAX_N = 1000;
  node_t V[MAX_N]; int N;
};

void Prim(const graph_t& graph, vector<int>& cost, vector<const edge_t*>& parent, int root) {
  cost.assign(graph.N, INFINITE); parent.assign(graph.N, (const edge_t*)NULL);
  vector<bool> mark(graph.N, false);
  typedef pair<int, int> state_t;  // (cost, node)
  priority_queue<state_t, vector<state_t>, greater<state_t> > pqueue;
  for (pqueue.push(state_t(cost[root] = 0, root)); !pqueue.empty();) {
    const int node = pqueue.top().second; pqueue.pop();
    if (mark[node]) continue; mark[node] = true;
    for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
      const edge_t* const edge = &graph.V[node][k];
      if (mark[edge->dest] || edge->cost >= cost[edge->dest]) continue;
      cost[edge->dest] = edge->cost, parent[edge->dest] = edge;
      pqueue.push(state_t(edge->cost, edge->dest));
    }
  }
}

};

////////////////////// disjktra ////////////////////////////

namespace disjktra {

const int INFINITE = 100000;
struct edge_t {
  edge_t(int src_, int dest_, int cost_)
    : src(src_), dest(dest_), cost(cost_) {}
  int src, dest, cost;
};
typedef vector<edge_t> node_t;
struct graph_t {
  static const int MAX_N = 500;
  node_t V[MAX_N]; int N;
};
struct path_t {
  vector<const edge_t*> edges;
  int cost;
};

void GetTreePath(const vector<const edge_t*>& parent, int dest, path_t& path) {
  if (parent[dest] != NULL) {
    GetTreePath(parent, parent[dest]->src, path);
    path.edges.push_back(parent[dest]);
    path.cost += parent[dest]->cost;
  }
}

void Disjktra(const graph_t& graph, vector<int>& cost, vector<const edge_t*>& parent, int root) {
  cost.assign(graph.N, INFINITE); parent.assign(graph.N, (const edge_t*)NULL);
  typedef pair<int, int> state_t;  // (cost, node)
  priority_queue<state_t, vector<state_t>, greater<state_t> > pqueue;
  for (pqueue.push(state_t(cost[root] = 0, root)); !pqueue.empty();) {
    const state_t state = pqueue.top(); pqueue.pop();
    const int nextCost = state.first, node = state.second;
    if (nextCost > cost[node]) continue;
    for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
      const edge_t* const edge = &graph.V[node][k];
      const int newCost = nextCost + edge->cost;
      if (newCost >= cost[edge->dest]) continue;
      cost[edge->dest] = newCost, parent[edge->dest] = edge;
      pqueue.push(state_t(newCost, edge->dest));
    }
  }
}

};

////////////////////// bellmann ford ////////////////////////////

namespace bellmann_ford {

const int INFINITE = 1000000;
struct edge_t {
  edge_t(int src_, int dest_, int cost_)
    : src(src_), dest(dest_), cost(cost_) {}
  int src, dest, cost;
};
typedef vector<edge_t> node_t;
struct graph_t {
  static const int MAX_N = 500;
  node_t V[MAX_N]; int N;
};
struct path_t {
  vector<const edge_t*> edges;
  int cost;
};

pair<path_t, bool> GetNegativeCycle(const vector<const edge_t*>& tree, int node) {
  vector<bool> mark(tree.size(), false); mark[node] = true;
  path_t path; path.cost = 0;
  while (true) {
    if (tree[node] == NULL) return make_pair(path, false);
    path.edges.push_back(tree[node]); path.cost += tree[node]->cost;
    if (mark[node = tree[node]->src]) break; mark[node] = true;
  }
  return make_pair(path, true);
}

pair<path_t, bool> BellmannFord(const graph_t& graph, vector<int>& cost,
                                vector<const edge_t*>& parent, int root) {
  pair<path_t, bool> cycle; cycle.second = false;
  queue<int> active; active.push(root); active.push(graph.N);
  vector<bool> mark(graph.N, false); mark[root] = true;
  cost.assign(graph.N, INFINITE); cost[root] = 0;
  parent.assign(graph.N, (const edge_t*)NULL);
  for (int pass = 0; pass < graph.N && !active.empty() && !cycle.second;) {
    const int node = active.front(); active.pop();
    if (node == graph.N && ++pass < graph.N) {
      active.push(graph.N);
    } else if (node < graph.N) {
      mark[node] = false;
      for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
        const edge_t* const edge = &graph.V[node][k];
        const int newCost = cost[node] + edge->cost;
        if (newCost < cost[edge->dest]) {
          cost[edge->dest] = newCost, parent[edge->dest] = edge;
          cycle = GetNegativeCycle(parent, edge->dest);
          if (!mark[edge->dest]) active.push(edge->dest), mark[edge->dest] = true;
        }
      }
    }
  }
  return cycle;
}

};

////////////////////// floyd warshall ////////////////////////////

namespace floyd_warshall {  // <-- pass to static graph

typedef vector<vector<int> > graph_t;
void FloydWarshall(const graph_t& graph, graph_t& distance, graph_t& parent) {
  int i, j, k, N = graph.size();
  for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    distance[i][j] = i == j ? 0 : graph[i][j], parent[i][j] = i;
  for (k = 0; k < N; ++k) for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    if (distance[i][k] >= 0 && distance[k][j] >= 0) {
      const int newPath = distance[i][k] + distance[k][j];
      if (distance[i][j] < 0 || newPath < distance[i][j])
        distance[i][j] = newPath, parent[i][j] = parent[k][j];
    }
}

void TransitiveClosure(const graph_t& graph, graph_t& hull, graph_t& parent) {
  int i, j, k, N = graph.size();
  for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    hull[i][j] = i == j || graph[i][j] != 0, parent[i][j] = i;
  for (k = 0; k < N; ++k) for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    if (!hull[i][j] && hull[i][k] && hull[k][j])
      hull[i][j] = true, parent[i][j] = parent[k][j];
}

void MiniMax(const graph_t& graph, graph_t& minimax, graph_t& parent) {
  int i, j, k, N = graph.size();
  for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    minimax[i][j] = i == j ? 0 : graph[i][j], parent[i][j] = i;
  for (k = 0; k < N; ++k) for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    if (minimax[i][k] >= 0 && minimax[k][j] >= 0) {
      const int newPath = max(minimax[i][k], minimax[k][j]);
      if (minimax[i][j] < 0 || newPath < minimax[i][j])
        minimax[i][j] = newPath, parent[i][j] = parent[k][j];
    }
}

void MaxiMin(const graph_t& graph, graph_t& maximin, graph_t& parent) {
  int i, j, k, N = graph.size();
  for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    maximin[i][j] = i == j ? 0 : graph[i][j], parent[i][j] = i;
  for (k = 0; k < N; ++k) for (i = 0; i < N; ++i) for (j = 0; j < N; ++j)
    if (maximin[i][k] >= 0 && maximin[k][j] >= 0) {
      const int newPath = min(maximin[i][k], maximin[k][j]);
      if (maximin[i][j] < 0 || newPath > maximin[i][j])
        maximin[i][j] = newPath, parent[i][j] = parent[k][j];
    }
}

void GetFloydPath(const graph_t& parent, int src, int dest, vector<int>& path) {
  if (src != dest) GetFloydPath(parent, src, parent[src][dest], path);
  path.push_back(dest);
}

};

////////////////////// preflow-push max-flow ////////////////////////////

namespace preflow_push {

const int INFINITE = 100000;
struct edge_t {
  edge_t(int src_, int dest_, int rev_, int capacity_) :
    src(src_), dest(dest_), rev(rev_), capacity(capacity_), flow(0) {}
  int src, dest, rev, capacity, flow;
};
struct node_t {
  vector<edge_t> edges; int excess;
  node_t() : excess(0) {}
};
struct graph_t {
  static const int MAX_N = 500;
  node_t V[MAX_N]; int N;
};
struct flow_t {
  flow_t(int N) : height(N, N), count(2 * N, 0) {count[N] = N;}
  vector<int> height, count;
};

void AddDirectedEdge(graph_t& graph, int src, int dest, int capacity) {
  const edge_t edge1(src, dest, graph.V[dest].edges.size(), capacity);
  const edge_t edge2(dest, src, graph.V[src].edges.size(), 0);
  graph.V[src].edges.push_back(edge1); graph.V[dest].edges.push_back(edge2);
}

int ChangeHeight(flow_t& data, int node, int newHeight) {
  const int oldHeight = data.height[node];
  ++data.count[data.height[node] = newHeight];
  return --data.count[oldHeight] <= 0 ? oldHeight : -1;
}

void InitHeights(const graph_t& graph, flow_t& data, int src, int dest) {
  queue<int> nqueue; nqueue.push(dest); ChangeHeight(data, dest, 0);
  while (!nqueue.empty()) {
    const int node = nqueue.front(); nqueue.pop();
    for (int k = 0, N = (int)graph.V[node].edges.size(); k < N; ++k) {
      const edge_t* const edge = &graph.V[node].edges[k];
      if (edge->dest != src && data.height[node] + 1 < data.height[edge->dest]) {
        ChangeHeight(data, edge->dest, data.height[node] + 1);
        nqueue.push(edge->dest);
      }
    }
  }
}

void AddFlow(graph_t& graph, edge_t* edge, int flow) {
  edge->capacity -= flow, edge->flow += flow;
  graph.V[edge->dest].excess += flow;
}

void PushFlow(graph_t& graph, edge_t* edge, int flow) {
  edge_t* const rev = &graph.V[edge->dest].edges[edge->rev];
  AddFlow(graph, edge, +flow); AddFlow(graph, rev, -flow);
}

int PreflowPushMaxFlow(graph_t& graph, int src, int dest) {
  flow_t data(graph.N); InitHeights(graph, data, src, dest);
  vector<bool> mark(graph.N, false); queue<int> active;
  for (int k = 0, N = (int)graph.V[src].edges.size(); k < N; ++k) {
    edge_t* const edge = &graph.V[src].edges[k];
    PushFlow(graph, edge, edge->capacity);
    if (edge->dest != src && edge->dest != dest && !mark[edge->dest])
      active.push(edge->dest), mark[edge->dest] = true;
  }
  while (!active.empty()) {
    const int node = active.front(); int minHeight = INFINITE;
    for (int k = 0, N = (int)graph.V[node].edges.size(); k < N && graph.V[node].excess > 0; ++k) {
      edge_t* const edge = &graph.V[node].edges[k]; if (edge->capacity <= 0) continue;
      if (data.height[node] > data.height[edge->dest]) {
        PushFlow(graph, edge, min(graph.V[node].excess, edge->capacity));
        if (edge->dest != src && edge->dest != dest && !mark[edge->dest])
          active.push(edge->dest), mark[edge->dest] = true;
      } else minHeight = min(minHeight, data.height[edge->dest]);
    }
    if (graph.V[node].excess > 0) {
      const int oldHeight = ChangeHeight(data, node, minHeight + 1);
      if (oldHeight != -1) {
        if (oldHeight > data.height[src]) break;
        for (int k = 0; k < graph.N; ++k) if (k != src && data.height[k] > oldHeight)
          ChangeHeight(data, node, max(data.height[k], data.height[src] + 1));
      }
    } else active.pop(), mark[node] = false;
  }
  return graph.V[dest].excess;
}

};

////////////////////// augmenting shortest-path max-flow ////////////////////////////

namespace augmenting_max_flow {

const int INFINITE = 100000;
struct edge_t {
  edge_t(int src_, int dest_, int rev_, int capacity_) :
    src(src_), dest(dest_), rev(rev_), capacity(capacity_), flow(0) {}
  int src, dest, rev, capacity, flow;
};
typedef vector<edge_t> node_t;
struct graph_t {
  static const int MAX_N = 500;
  node_t V[MAX_N]; int N;
};
struct path_t {
  vector<edge_t*> edges;
  vector<int> capacities;
  int flow;
};
struct flow_t {
  flow_t(int N) : distance(N, N), next(N, 0), count(N + 1, 0) {count[N] = N;}
  vector<int> distance, next, count;
};

void AddDirectedEdge(graph_t& graph, int src, int dest, int capacity) {
  const edge_t edge1(src, dest, (int)graph.V[dest].size(), capacity);
  const edge_t edge2(dest, src, (int)graph.V[src].size(), 0);
  graph.V[src].push_back(edge1); graph.V[dest].push_back(edge2);
}

void ReverseBreadthFirstSearch(graph_t& graph, flow_t& data, int src) {
  data.distance[src] = 0, ++data.count[0], --data.count[graph.N];
  queue<int> nqueue; nqueue.push(src);
  while (!nqueue.empty()) {
    const int node = nqueue.front(); nqueue.pop();
    for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
      if (data.distance[graph.V[node][k].dest] < graph.N) continue;
      --data.count[data.distance[graph.V[node][k].dest]];
      data.distance[graph.V[node][k].dest] = data.distance[node] + 1;
      ++data.count[data.distance[graph.V[node][k].dest]];
      nqueue.push(graph.V[node][k].dest);
    }
  }
}

void AddFlow(edge_t& edge, int flow) {
  edge.capacity -= flow;
  edge.flow += flow;
}

void Augment(graph_t& graph, path_t& path) {
  for (int k = 0, N = (int)path.edges.size(); k < N; k++) {
    edge_t& edge = *path.edges[k]; AddFlow(edge, +path.flow);
    AddFlow(graph.V[edge.dest][edge.rev], -path.flow);
  }
  path.edges.clear(); path.capacities.clear();
  path.flow = INFINITE;
}

bool Advance(graph_t& graph, flow_t& data, path_t& path, int& node) {
  bool found = false;
  for (int k = data.next[node], N = (int)graph.V[node].size(); k < N && !found; ++k) {
    if (graph.V[node][k].capacity <= 0) continue;
    if (data.distance[node] == data.distance[graph.V[node][k].dest] + 1) {
      path.edges.push_back(&graph.V[node][k]);
      path.capacities.push_back(path.flow);
      path.flow = min(path.flow, graph.V[node][k].capacity);
      data.next[node] = k, node = graph.V[node][k].dest, found = true;
    }
  }
  if (!found) data.next[node] = 0; return found;
}

bool Retreat(graph_t& graph, flow_t& data, path_t& path, int& node) {
  int closestNode = graph.N - 1;
  for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
    if (graph.V[node][k].capacity <= 0) continue;
    closestNode = min(closestNode, data.distance[graph.V[node][k].dest]);
  }
  if (--data.count[data.distance[node]] <= 0) return false;
  data.distance[node] = closestNode + 1, ++data.count[data.distance[node]];
  if (!path.edges.empty()) {
    node = path.edges.back()->src; path.edges.pop_back();
    path.flow = path.capacities.back(); path.capacities.pop_back();
  }
  return true;
}

int AugmentingPathMaxFlow(graph_t& graph, int src, int dest) {
  flow_t data(graph.N);
  ReverseBreadthFirstSearch(graph, data, dest);
  path_t path; path.flow = INFINITE;
  int networkFlow = 0, node = src;
  while (data.distance[src] < graph.N) {
    if (node == dest) {
      networkFlow += path.flow; Augment(graph, path); node = src;
    } else if (!Advance(graph, data, path, node)) {
      if (!Retreat(graph, data, path, node)) break;
    }
  }
  return networkFlow;
}

};

////////////////////// successive shortest path min-cost max-flow ////////////////////////////

namespace min_cost_max_flow {

const int INFINITE = 100000;
struct edge_t {
  edge_t(int src_, int dest_, int rev_, int capacity_, int cost_) :
    src(src_), dest(dest_), rev(rev_), capacity(capacity_), cost(cost_), flow(0) {}
  int src, dest, rev, capacity, cost, flow;
};
typedef vector<edge_t> node_t;
struct graph_t {
  static const int MAX_N = 700;
  node_t V[MAX_N]; int N;
};
struct path_t {
  vector<edge_t*> edges;
  int flow, cost;
};

void AddDirectedEdge(graph_t& graph, int src, int dest, int capacity, int cost) {
  const edge_t edge1(src, dest, (int)graph.V[dest].size(), capacity, cost);
  const edge_t edge2(dest, src, (int)graph.V[src].size(), 0, -cost);
  graph.V[src].push_back(edge1); graph.V[dest].push_back(edge2);
}

void BellmannFord(graph_t& graph, vector<edge_t*>& tree, vector<int>& distance, int src) {
  queue<int> active; active.push(src); active.push(graph.N);
  vector<bool> mark(graph.N, false); mark[src] = true;
  tree[src] = NULL; distance[src] = 0;
  for (int pass = 0; pass < graph.N && !active.empty();) {
    const int node = active.front(); active.pop();
    if (node == graph.N && ++pass < graph.N) {
      active.push(graph.N);
    } else if (node < graph.N) {
      mark[node] = false;
      for (int j = 0, N = (int)graph.V[node].size(); j < N; ++j) {
        edge_t* const edge = &graph.V[node][j]; if (edge->capacity <= 0) continue;
        const int newDistance = distance[node] + edge->cost;
        if (newDistance < distance[edge->dest]) {
          tree[edge->dest] = edge, distance[edge->dest] = newDistance;
          if (!mark[edge->dest]) active.push(edge->dest), mark[edge->dest] = true;
        }
      }
    }
  }
}

void Disjktra(graph_t& graph, vector<edge_t*>& tree, vector<int>& distance,
              const vector<int>& potential, int src) {
  typedef pair<int, pair<int, edge_t*> > state_t;  // (distance, (node, parent))
  priority_queue<state_t, vector<state_t>, greater<state_t> > pqueue;
  pqueue.push(state_t(0, make_pair(src, (edge_t*)NULL)));
  while (!pqueue.empty()) {
    const int nextDistance = pqueue.top().first, node = pqueue.top().second.first;
    edge_t* const pred = pqueue.top().second.second; pqueue.pop();
    if (nextDistance >= distance[node]) continue;
    distance[node] = nextDistance, tree[node] = pred;
    for (int k = 0, N = (int)graph.V[node].size(); k < N; ++k) {
      edge_t* const edge = &graph.V[node][k]; if (edge->capacity <= 0) continue;
      const int newDistance = nextDistance + edge->cost + potential[edge->dest] - potential[node];
      if (newDistance < distance[edge->dest])
        pqueue.push(state_t(newDistance, make_pair(edge->dest, edge)));
    }
  }
}

bool GetPath(const vector<edge_t*>& parent, path_t& path, int src, int dest) {
  path.edges.clear(); path.flow = INFINITE, path.cost = 0;
  while (src != dest) {
    if (parent[dest] == NULL) return false;
    path.edges.push_back(parent[dest]);
    path.flow = min(path.flow, parent[dest]->capacity);
    path.cost += parent[dest]->cost;
    dest = parent[dest]->src;
  }
  return true;
}

void AddFlow(edge_t& edge, int flow) {
  edge.capacity -= flow;
  edge.flow += flow;
}

void Augment(graph_t& graph, path_t& path) {
  for (int k = 0, N = (int)path.edges.size(); k < N; k++) {
    edge_t& edge = *path.edges[k]; AddFlow(edge, +path.flow);
    AddFlow(graph.V[edge.dest][edge.rev], -path.flow);
  }
}

pair<int, int> MinCostMaxFlow(graph_t& graph, int src, int dest) {
  vector<edge_t*> tree(graph.N, (edge_t*)NULL);
  vector<int> potential(graph.N, 0), distance(graph.N, INFINITE);
  BellmannFord(graph, tree, distance, src);
  for (int k = 0; k < graph.N; ++k) potential[k] = -distance[k];
  int networkFlow = 0, flowCost = 0;
  for (path_t path; GetPath(tree, path, src, dest);) {
    networkFlow += path.flow, flowCost += path.cost * path.flow; Augment(graph, path);
    for (int k = 0; k < graph.N; ++k) tree[k] = NULL, distance[k] = INFINITE;
    Disjktra(graph, tree, distance, potential, src);
    for (int k = 0; k < graph.N; ++k) potential[k] -= distance[k];
  }
  return make_pair(networkFlow, flowCost);
}

};

////////////////////// longest increasing sequence ////////////////////////////

namespace lis {

vector<int> GetSequence(const vector<int>& predecesor, int index) {
  vector<int> sequence;
  while (index != -1) {
    sequence.push_back(index);
    index = predecesor[index];
  }
  reverse(sequence.begin(), sequence.end());
  return sequence;
}

namespace ON2 {

template <class Iterator, class Pred>
vector<int> LongestSequence(Iterator begin, Iterator end, Pred pred) {  // O(N^2)
  vector<int> length(end - begin, 1), predecesors(end - begin, -1);
  int best = 0;
  for (Iterator next = begin + 1; next != end; ++next) {
    for (Iterator prev = begin; prev != next; ++prev)
      if (pred(*prev, *next) && length[prev - begin] + 1 > length[next - begin]) {
        length[next - begin] = length[prev - begin] + 1;
        predecesors[next - begin] = prev - begin;
      }
    if (length[next - begin] > length[best])
      best = next - begin;
  }
  return GetSequence(predecesors, best);
}

};

namespace ONlgK {

template <class Iterator, class Pred>
struct IndexFunctor {
  Iterator begin; const Pred& pred;
  IndexFunctor(Iterator begin_, const Pred& pred_) : begin(begin_), pred(pred_) {}
  bool operator ()(int index1, int index2) {
    return pred(*(begin + index1), *(begin + index2));
  }
};

template <class Iterator, class Pred>
vector<int> LongestSequence(Iterator begin, Iterator end, Pred pred) {  // O(N*lgK)
  vector<int> elements(1, 0), length(end - begin, 1), predecesors(end - begin, -1);
  int best = 0;
  for (Iterator next = begin + 1; next != end; ++next) {
    const int index = next - begin;
    vector<int>::iterator ptr = lower_bound(elements.begin(), elements.end(), index,
                                            IndexFunctor(begin, pred));
    if (ptr != elements.end()) *ptr = index;
    else elements.push_back(index), ptr = elements.end() - 1;
    if (ptr != elements.begin()) {
      const int pred = predecesors[index] = ptr[-1];
      if ((length[index] = length[pred] + 1) > length[best]) best = index;
    }
  }
  return GetSequence(predecesors, best);
}

};

};

////////////////////// longest common subsequence ////////////////////////////

namespace lcs {

enum {ROW = 0, COL, BOTH};
vector<int> GetSequence(const vector<vector<int> >& predecesor, int i, int j) {
  vector<int> sequence;
  while (i > 0 && j > 0) {
    if (predecesor[i][j] == ROW) --i;
    else if (predecesor[i][j] == COL) --j;
    else sequence.push_back(--i), --j;
  }
  reverse(all(sequence));
  return sequence;
}

template <class Iterator, class Pred>
vector<int> LongestCommonSubsequence(Iterator begin1, Iterator end1, Iterator begin2, Iterator end2, Pred pred) {
  const int length1 = end1 - begin1, length2 = end2 - begin2;
  vector<vector<int> > length(length1 + 1, vector<int>(length2 + 1, 0));
  vector<vector<int> > predecesor(length1 + 1, vector<int>(length2 + 1, -1));
  for (int i = 1; i <= length1; ++i) for (int j = 1; j <= length2; ++j) {
    if (pred(begin1[i - 1], begin2[j - 1]))
      length[i][j] = length[i - 1][j - 1] + 1, predecesor[i][j] = BOTH;
    else if (length[i - 1][j] >= length[i][j - 1])
      length[i][j] = length[i - 1][j], predecesor[i][j] = ROW;
    else length[i][j] = length[i][j - 1], predecesor[i][j] = COL;
  }
  return GetSequence(predecesor, end1 - begin1, end2 - begin2);
}

};

////////////////////// minimum edit distance ////////////////////////////

enum {END = -1, EQUAL = 0, MODIFY, INSERT, DELETE};
typedef pair<int, pair<int, int> > edition_t;  // (action, (src, index))
int GetEditions(const vector<vector<int> >& previous, int i, int j, vector<edition_t>& editions) {
  if (previous[i][j] == EQUAL) return GetEditions(previous, i - 1, j - 1, editions);
  else if (previous[i][j] == MODIFY) {
    const int count = GetEditions(previous, i - 1, j - 1, editions);
    editions.push_back(edition_t(MODIFY, make_pair(j - 1, j - 1)));
    return count;
  } else if (previous[i][j] == INSERT) {
    const int count = GetEditions(previous, i, j - 1, editions);
    editions.push_back(edition_t(INSERT, make_pair(j - 1, i - count)));
    return count - 1;
  } else if (previous[i][j] == DELETE) {
    const int count = GetEditions(previous, i - 1, j, editions);
    editions.push_back(edition_t(DELETE, make_pair(i - 1, i - count - 1)));
    return count + 1;
  } else return 0;
}

template <class T, class Pred>
vector<edition_t> MinimumEditDistance(const vector<T>& V1, const vector<T>& V2, Pred pred) {
  vector<vector<int> > distance(V1.size() + 1, vector<int>(V2.size() + 1));
  vector<vector<int> > previous(V1.size() + 1, vector<int>(V2.size() + 1));
  distance[0][0] = 0, previous[0][0] = END;
  for (int i = 1; i <= V1.size(); ++i) distance[i][0] = i, previous[i][0] = DELETE;
  for (int j = 1; j <= V2.size(); ++j) distance[0][j] = j, previous[0][j] = INSERT;
  for (int i = 1; i <= V1.size(); ++i) for (int j = 1; j <= V2.size(); ++j) {
    if (pred(V1[i - 1], V2[j - 1])) {
      distance[i][j] = distance[i - 1][j - 1], previous[i][j] = EQUAL;
    } else {
      const int modDist = distance[i - 1][j - 1] + 1;
      const int insDist = distance[i][j - 1] + 1;
      const int delDist = distance[i - 1][j] + 1;
      if (modDist <= insDist && modDist <= delDist)
        distance[i][j] = modDist, previous[i][j] = MODIFY;
      else if (insDist <= modDist && insDist <= delDist)
        distance[i][j] = insDist, previous[i][j] = INSERT;
      else distance[i][j] = delDist, previous[i][j] = DELETE;
    }
  }
  vector<edition_t> editions;
  GetEditions(previous, V1.size(), V2.size(), editions);
  return editions;
}

////////////////////// zero-one knapsack ////////////////////////////

namespace _0_1_knapsack {

vector<int> GetKnapsack(const vector<vector<int> >& predecesor, int i, int w) {
  vector<int> elements;
  for (; i > 0; --i) if (predecesor[i][w] != -1)
    elements.push_back(i - 1), w -= predecesor[i][w];
  reverse(elements.begin(), elements.end());
  return elements;
}

vector<int> ZeroOneKnapsack(const vector<int>& price, const vector<int>& weight, int maxWeight) {
  vector<vector<int> > value(price.size() + 1, vector<int>(maxWeight + 1, 0));
  vector<vector<int> > predecesor(price.size() + 1, vector<int>(maxWeight + 1, -1));
  for (int i = 1; i <= price.size(); ++i) {
    for (int w = 0; w < weight[i - 1] && w <= maxWeight; ++w)
      value[i][w] = value[i - 1][w];
    for (int w = weight[i - 1]; w <= maxWeight; ++w) {
      const int newPrice = price[i - 1] + value[i - 1][w - weight[i - 1]];
      if (newPrice > value[i - 1][w])
        value[i][w] = newPrice, predecesor[i][w] = weight[i - 1];
      else value[i][w] = value[i - 1][w];
    }
  }
  return GetKnapsack(predecesor, price.size(), maxWeight);
}

};

////////////////////// coin change ////////////////////////////

vector<int> InfiniteChangeMinCoin(const vector<int>& coin, const int maxMoney) {
  vector<int> minCoin(maxMoney + 1, -1); minCoin[0] = 0;
  for (int k = 0; k < coin.size(); k++) for (int i = coin[k]; i <= maxMoney; i++) {
    if (minCoin[i - coin[k]] == -1) continue;
    const int neededCoins = minCoin[i - coin[k]] + 1;
    if (minCoin[i] == -1 || neededCoins < minCoin[i])
      minCoin[i] = neededCoins;
  }
  return minCoin;
}

vector<int> LimitedChangeMinCoin(const vector<int>& coin, const int maxMoney) {
  vector<int> minCoin(maxMoney + 1, -1); minCoin[0] = 0;
  vector<int> tmpCoin(maxMoney + 1, -1); tmpCoin[0] = 0;
  for (int k = 0; k < coin.size(); ++k) {
    for (int i = coin[k]; i <= maxMoney; ++i) if (minCoin[i - coin[k]] != -1) {
      const int neededCoins = minCoin[i - coin[k]] + 1;
      if (tmpCoin[i] == -1 || neededCoins < tmpCoin[i])
        tmpCoin[i] = neededCoins;
    }
    copy(tmpCoin.begin(), tmpCoin.end(), minCoin.begin());
  }
  return minCoin;
}

////////////////////// optimal array multiplication ////////////////////////////

typedef vector<vector<int> > matrix_t;
pair<int, int> ProcessTree(const vector<int>& A, const matrix_t& pred, int i, int j) {
  if (i == j) return make_pair(A[i], A[i + 1]);
  else {
    const pair<int, int> left = ProcessTree(A, pred, i, pred[i][j]);
    const pair<int, int> right = ProcessTree(A, pred, pred[i][j] + 1, j);
    assert(left.second == right.first);
    return make_pair(left.first, right.second);
  }
}

void OptimalArrayMultiplication(const vector<int>& A, matrix_t& cost, matrix_t& pred) {
  const int N = A.size() - 1;
  cost = matrix_t(N, vector<int>(N, 0));
  pred = matrix_t(N, vector<int>(N, -1));
  for (int k = 1; k < N; ++k) for (int i = 0; i < N - k; ++i) {
    const int j = i + k;
    cost[i][j] = cost[i][i] + cost[i + 1][j] + A[i] * A[i + 1] * A[j + 1];
    pred[i][j] = i;
    for (int m = i + 1; m < j; ++m) {
      const int newCost = cost[i][m] + cost[m + 1][j] + A[i] * A[m + 1] * A[j + 1];
      if (newCost < cost[i][j]) cost[i][j] = newCost, pred[i][j] = m;
    }
  }
}

////////////////////// orientation and triangles ////////////////////////////

namespace new_geometry {  // erase this namespace once all the geometry stuff is ported to templates

enum {NONE = 0, POINT, LINE, TANGENT, SECANT, CIRCLE};
template <class T> struct point_t {
  T x, y; point_t() {}
  point_t(T x_, T y_) : x(x_), y(y_) {}
};

template <class T>
T TriangleArea(const point_t<T>& a, const point_t<T>& b, const point_t<T>& c) {
  // calculates the double of the signed area defined by the points ABC
  return ((c.y - a.y) * (b.x - a.x) - (b.y - a.y) * (c.x - a.x));
}

// CCW --> +1, CW --> -1, Collinear --> 0
template <class T>
int Orientation(const point_t<T>& a, const point_t<T>& b, const point_t<T>& c, bool patch = false) {
  const T dx1 = b.x - a.x, dy1 = b.y - a.y;
  const T dx2 = c.x - a.x, dy2 = c.y - a.y;
  if (patch) {
    if (dx1 * dx2 < 0 || dy1 * dy2 < 0) return +1;
    if (dx1 * dx1 + dy1 * dy1 > dx2 * dx2 + dy2 * dy2) return -1;
  }
  if (dy1 * dx2 < dx1 * dy2) return +1;
  if (dy1 * dx2 > dx1 * dy2) return -1;
  return 0;
}

template <class T>
double InverseInterpolation(const point_t<T>& p0, const point_t<T>& p1, const point_t<T>& p2) {
  if (abs(p0.x - p1.x) > abs(p0.y - p1.y))
    return (p2.x - p0.x) / (double)(p1.x - p0.x);
  else return (p2.y - p0.y) / (double)(p1.y - p0.y);
}

template <class T>
bool IsInsideTriangle(const point_t<T>& point, const point_t<T>* triangle) {
  const T areaAB = abs(TriangleArea(triangle[0], point, triangle[1]));
  const T areaBC = abs(TriangleArea(triangle[1], point, triangle[2]));
  const T areaCA = abs(TriangleArea(triangle[2], point, triangle[0]));
  const T areaT = abs(TriangleArea(triangle[0], triangle[1], triangle[2]));
  return areaAB + areaBC + areaCA == areaT;
}

////////////////////// spheric distance ////////////////////////////

double haversine(double x) {
  return (1 - cos(x)) / 2;
//  return sin(x / 2) * sin(x / 2);
}

double SphericDistance(double lat1, double lon1, double lat2, double lon2, double radius) {
  const double a = haversine(lat2 - lat1);
  const double b = cos(lat1) * cos(lat2) * haversine(lon2 - lon1);
  const double c = 2 * atan2(sqrt(a + b + eps), sqrt(1 - a - b + eps));
  return radius * c;
}

////////////////////// line and segment intersections ////////////////////////////

template <class T>
bool CollinearIntersects(T x11, T x21, T x12, T x22) {
  if (x11 > x21) swap(x11, x21); if (x12 > x22) swap(x12, x22);
  if (x11 > x12) swap(x11, x12), swap(x21, x22);
  return x21 >= x12;  /* equal result */
}

template <class T> struct segment_t {
  point_t<T> p0, p1;
  segment_t(point_t<T> p0_, point_t<T> p1_) : p0(p0_), p1(p1_) {}
};

template <class T>
bool PointInSegment(const segment_t<T>& line, const point_t<T>& point) {
  if (!(line.p0.x <= point.x && point.x <= line.p1.x) &&
      !(line.p1.x <= point.x && point.x <= line.p0.x)) return false;
  if (!(line.p0.y <= point.y && point.y <= line.p1.y) &&
      !(line.p1.y <= point.y && point.y <= line.p0.y)) return false;
  return (line.p1.y - line.p0.y) * point.x - line.p1.y * line.p0.x ==
         (line.p1.x - line.p0.x) * point.y - line.p1.x * line.p0.y;
}

template <class T>
T PointToSegmentDistance(const point_t<T>& p, const segment_t<T>& s) {
  // return the square of the distance of the closest point in s to p
  const T dx = s.p1.x - s.p0.x, dy = s.p1.y - s.p0.y;
  const T dx1 = p.x - s.p0.x, dy1 = p.y - s.p0.y;
  const T dx2 = p.x - s.p1.x, dy2 = p.y - s.p1.y;
  const T dot1 = +dx * dx1 + dy * dy1, dot2 = -dx * dx2 - dx * dy2;
  if (dot1 < T(0)) return dx1 * dx1 + dy1 * dy1;
  if (dot2 < T(0)) return dx2 * dx2 + dy2 * dy2;
  const T dist2 = dx * dx + dy * dy;
  const T rx = s.p0.x + (dx * dot1) / dist2;
  const T ry = s.p0.y + (dy * dot1) / dist2;
  const T drx = rx - p.x, dry = ry - p.y;
  return drx * drx + dry * dry;
}

template <class T>
bool Intersects(const segment_t<T>& seg1, const segment_t<T>& seg2) {
  int ccw1 = Orientation(seg2.p0, seg1.p0, seg2.p1, false);
  int ccw2 = Orientation(seg2.p0, seg1.p1, seg2.p1, false);
  int ccw3 = Orientation(seg1.p0, seg2.p0, seg1.p1, false);
  int ccw4 = Orientation(seg1.p0, seg2.p1, seg1.p1, false);
  if (ccw1 != 0 || ccw2 != 0 || ccw3 != 0 || ccw4 != 0)
    return ccw1 * ccw2 <= 0 && ccw3 * ccw4 <= 0;
  ccw1 = Orientation(seg2.p0, seg1.p0, seg2.p1, true);
  ccw2 = Orientation(seg2.p0, seg1.p1, seg2.p1, true);
  ccw3 = Orientation(seg1.p0, seg2.p0, seg1.p1, true);
  ccw4 = Orientation(seg1.p0, seg2.p1, seg1.p1, true);
  return ccw1 * ccw2 <= 0 && ccw3 * ccw4 <= 0;
}

template <class T> struct line_t {
  point_t<T> p0, p1;
  line_t(point_t<T> p0_, point_t<T> p1_) : p0(p0_), p1(p1_) {}
};
template <class T> struct plane_t {
  T a, b, c;
  plane_t() : a(0), b(0), c(0) {}
  plane_t(T a_, T b_, T c_) : a(a_), b(b_), c(c_) {}
};

template <class T>
pair<int, point_t<T> > IntersectLines(const line_t<T>& line0, const line_t<T>& line1) {
  // calculates the intersection between two infinite lines
  const T dx0 = line0.p1.x - line0.p0.x, dx1 = line1.p1.x - line1.p0.x;
  const T dy0 = line0.p1.y - line0.p0.y, dy1 = line1.p1.y - line1.p0.y;
  const T numT0 = dy1 * (line0.p0.x - line1.p0.x) - dx1 * (line0.p0.y - line1.p0.y);
  const T denT0 = dx1 * dy0 - dy1 * dx0;
  if (denT0 == T(0)) {
    if (numT0 == T(0)) return pair<int, point_t<T> >(LINE, point_t<T>(T(0), T(0)));
    else return pair<int, point_t<T> >(NONE, point_t<T>(T(0), T(0)));
  }
  const T t0 = numT0 / denT0, x = line0.p0.x + dx0 * t0, y = line0.p0.y + dy0 * t0;
  return pair<int, point_t<T> >(POINT, point_t<T>(x, y));
}

template <class T>
pair<int, point_t<T> > LinePlaneIntersection(const plane_t<T>& plane, const line_t<T>& line) {
  // calculates the intersection between a plane and an infinite line
  const T dx = line.p1.x - line.p0.x, dy = line.p1.y - line.p0.y;
  const T numT = plane.a * line.p0.x + plane.b * line.p0.y + plane.c;
  const T denT = plane.a * -dx + plane.b * -dy;
  if (denT == T(0)) {
    if (numT == T(0)) return pair<int, point_t<T> >(LINE, point_t<T>(T(0), T(0)));
    else return pair<int, point_t<T> >(NONE, point_t<T>(T(0), T(0)));
  }
  const T t = numT / denT, x = line.p0.x + t * dx, y = line.p0.y + t * dy;
  return pair<int, point_t<T> >(POINT, point_t<T>(x, y));
}

////////////////////// circles ////////////////////////////

template <class T> struct circle_t {
  circle_t(const point_t<T>& c_, T r2_)
    : c(c_), r2(r2_) {}
  point_t<T> c; T r2;
};

template <class T>
pair<circle_t<T>, bool> GetCircle(const point_t<T>& p1, const point_t<T>& p2, const point_t<T>& p3) {
  const T dx2x1 = p2.x - p1.x, dx2x3 = p2.x - p3.x;
  const T dy2y1 = p2.y - p1.y, dy2y3 = p2.y - p3.y;
  const T halfDen = dy2y3 * dx2x1 - dy2y1 * dx2x3;
  if (halfDen == 0) return pair<circle_t<T>, bool>(circle_t<T>(point_t<T>(T(0), T(0)), 0), false);
  const T d2x1 = p1.x * p1.x, d2x2 = p2.x * p2.x, d2x3 = p3.x * p3.x;
  const T d2y1 = p1.y * p1.y, d2y2 = p2.y * p2.y, d2y3 = p3.y * p3.y;
  const T coef1 = (d2x1 - d2x2) + (d2y1 - d2y2), coef2 = (d2x3 - d2x2) + (d2y3 - d2y2);
  const T cx = (coef2 * dy2y1 - coef1 * dy2y3) / (halfDen * T(2));
  const T cy = (coef1 * dx2x3 - coef2 * dx2x1) / (halfDen * T(2));
  const T dx1cx = p1.x - cx, dy1cy = p1.y - cy;
  const T r2 = dx1cx * dx1cx + dy1cy * dy1cy;
  return pair<circle_t<T>, bool>(circle_t<T>(point_t<T>(cx, cy), r2), true);
}

};

enum {NONE = 0, POINT, LINE, TANGENT, SECANT, CIRCLE};
typedef new_geometry::point_t<double> point_t;
typedef new_geometry::circle_t<double> circle_t;

double GetRadius(double a, double b, double c) {
  // calculate the radius that contains a triangle of sides {a, b, c}
  const double S = (a + b + c) / 2.0;
  const double L = sqrt(S * (S - a) * (S - b) * (S - c));
  return (a * b * c) / (4.0 * L);
}

int LineIntersectsCircle(const point_t& p1, const point_t& p2, const circle_t& c) {
  // check if an infinite line intersects with a circle
  const double x1 = p1.x - c.c.x, y1 = p1.y - c.c.y;
  const double x2 = p2.x - c.c.x, y2 = p2.y - c.c.y;
  const double dx = x2 - x1, dy = y2 - y1;
  const double dr2 = dx * dx + dy * dy;
  const double D = x1 * y2 - x2 * y1;
  const double delta2 = c.r2 * dr2 - D * D;
  if (DoubleComp(delta2, 0.0) == 0) return TANGENT;
  else if (DoubleComp(delta2, 0.0) < 0) return NONE;
  else return SECANT;
}

typedef pair<point_t, point_t> point_pair_t;
pair<int, point_pair_t> GetLineCircleIntersection(
    const point_t& p1, const point_t& p2, const circle_t& c) {
  // calculates the intersections between an infinite line and a circle
  const double x1 = p1.x - c.c.x, y1 = p1.y - c.c.y;
  const double x2 = p2.x - c.c.x, y2 = p2.y - c.c.y;
  const double dx = x2 - x1, dy = y2 - y1;
  const double dr2 = dx * dx + dy * dy;
  const double D = x1 * y2 - x2 * y1, delta2 = c.r2 * dr2 - D * D;
  if (DoubleComp(delta2, 0.0) < 0) return make_pair(NONE, point_pair_t());
  const double delta = sqrt(delta2 + eps);
  const double rx = DoubleSign(dy) * dx * delta, ry = fabs(dy) * delta;
  const double sx = D * dy, sy = -D * dx;
  point_pair_t intersection;
  intersection.first = point_t((sx - rx) / dr2 + c.c.x, (sy - ry) / dr2 + c.c.y);
  intersection.second = point_t((sx + rx) / dr2 + c.c.x, (sy + ry) / dr2 + c.c.y);
  return make_pair(DoubleComp(delta2, 0.0) == 0 ? TANGENT : SECANT, intersection);
}

pair<int, point_pair_t> GetTangentLines(const point_t& p, const circle_t& c) {
  // calculate tangent lines between a point and a circle
  const double dx = p.x - c.c.x, dy = p.y - c.c.y;
  const double dr2 = dx * dx + dy * dy, dr = sqrt(dr2 + eps), dt2 = dr2 - c.r2;
  if (DoubleComp(dt2, 0.0) < 0) return make_pair(0, point_pair_t());
  const double dt = sqrt(dt2 + eps);
  const double cosa = dt / dr, sina = sqrt(c.r2) / dr;
  const double vx = -dx * (dt / dr), vy = -dy * (dt / dr);
  point_t tangent1(cosa * vx - sina * vy,  sina * vx + cosa * vy);
  point_t tangent2(cosa * vx + sina * vy, -sina * vx + cosa * vy);
  tangent1.x += p.x, tangent1.y += p.y; tangent2.x += p.x, tangent2.y += p.y;
  return make_pair(DoubleComp(dt2, 0.0) == 0 ? 1 : 2, point_pair_t(tangent1, tangent2));
}

int CirclesIntersects(const circle_t& c1, const circle_t& c2) {
  // checks the number of intersections between both circles
  const double dx = c2.c.x - c1.c.x, dy = c2.c.y - c1.c.y;
  const double d = sqrt(dx * dx + dy * dy + eps);
  if (DoubleComp(d, 0.0) == 0 && DoubleComp(c1.r2, c2.r2) == 0) return CIRCLE;
  const double a = (c1.r2 - c2.r2 + d * d) / (2.0 * d);
  const double disc = c1.r2 - a * a;
  if (DoubleComp(disc, 0.0) == 0) return TANGENT;
  else if (DoubleComp(disc, 0.0) < 0) return NONE;
  else return SECANT;
}

pair<int, point_pair_t> GetCircleIntersection(const circle_t& c1, const circle_t& c2) {
  // calculate the intersection points between both circles
  const double dx = c2.c.x - c1.c.x, dy = c2.c.y - c1.c.y;
  const double d = sqrt(dx * dx + dy * dy + eps);
  if (DoubleComp(d, 0.0) == 0 && DoubleComp(c1.r2, c2.r2) == 0)
    return make_pair(CIRCLE, point_pair_t());
  const double a = (c1.r2 - c2.r2 + d * d) / (2.0 * d);
  const double disc = c1.r2 - a * a;
  if (DoubleComp(disc, 0.0) < 0) return make_pair(NONE, point_pair_t());
  const double h = sqrt(disc + eps), rx = -dy * (h / d), ry = dx * (h / d);
  const point_t center(c1.c.x + dx * (a / d), c1.c.y + dy * (a / d));
  const point_t point1(center.x - rx, center.y - ry), point2(center.x + rx, center.y + ry);
  return make_pair(DoubleComp(disc, 0.0) == 0 ? TANGENT : SECANT, point_pair_t(point1, point2));
}

////////////////////// point functors ////////////////////////////

namespace new_geometry {

template <class T>
struct PointSort : public binary_function<point_t<T>, point_t<T>, bool > {
  bool operator ()(const point_t<T>& p1, const point_t<T>& p2) {
    if (p1.x != p2.x) return p1.x < p2.x;
    else return p1.y < p2.y;
  }
};

template <class T>
struct PointEqual : public binary_function<point_t<T>, point_t<T>, bool > {
  bool operator ()(const point_t<T>& p1, const point_t<T>& p2) {
    return p1.x == p2.x && p1.y == p2.y;
  }
};

////////////////////// polygon area ////////////////////////////

template <class T> struct polygon_t {
  static const int MAX_N = 1000;
  polygon_t() : N(0) {}
  point_t<T> P[MAX_N]; int N;
};

template <class T>
T GetPolygonArea(const polygon_t<T>& polygon) {
  // calculates and returns the double of the signed area of the polygon
  T area = 0;
  for (int k = 0, N = polygon.N; k < N; ++k) {
    const point_t<T>& p0 = polygon.P[(k + 0) % N];
    const point_t<T>& p1 = polygon.P[(k + 1) % N];
    area += p0.x * p1.y - p1.x * p0.y;
  }
  return area;
}

////////////////////// center of mass ////////////////////////////

template <class T>
point_t<T> GetCenterOfMass(const polygon_t<T>& polygon) {
  const T area6 = GetPolygonArea(polygon) * 3; point_t<T> center(0, 0);
  for (int k = 0, N = polygon.N; k < N; ++k) {
    const point_t<T>& p1 = polygon[(k + 0) % N];
    const point_t<T>& p2 = polygon[(k + 1) % N];
    center.x += (p1.x + p2.x) * (p1.x * p2.y - p2.x * p1.y);
    center.y += (p1.y + p2.y) * (p1.x * p2.y - p2.x * p1.y);
  }
  center.x /= area6, center.y /= area6; return center;
}

};

////////////////////// closest pair ////////////////////////////

double EuclidDistance(const point_t& p1, const point_t& p2) {
  const double dx = p2.x - p1.x, dy = p2.y - p1.y;
  return sqrt(dx * dx + dy * dy);
}

pair<int, int> PlaneSweep(vector<point_t>& points) {
  sort(all(points), new_geometry::PointSort<double>());
  pair<int, int> result(0, 1);
  double minDist = EuclidDistance(points[0], points[1]);
  for (int k = 2; k < points.size(); ++k) {
    for (int i = k - 1; i >= 0; --i) {
      if (points[k].x - points[i].x > minDist) break;
      const double dist = EuclidDistance(points[k], points[i]);
      if (dist < minDist) minDist = dist, result = make_pair(i, k);
    }
  }
  return result;
}

////////////////////// pick's theorem ////////////////////////////

namespace new_geometry {

template <class T>
T GetBoundaryPoints(const polygon_t<T>& polygon) {
  T boundary = 0;
  for (int k = 0; k < polygon.N; ++k) {
    const point_t<T>& p1 = polygon.P[(k + 0) % polygon.N];
    const point_t<T>& p2 = polygon.P[(k + 1) % polygon.N];
    const T dx = p2.x - p1.x, dy = p2.y - p1.y;
    boundary += GCD(dx, dy);
  }
  return boundary;
}

template <class T>
i64 GetInteriorPoints(const polygon_t<T>& polygon) {
  const T area2 = GetPolygonArea(polygon);
  const T boundary = GetBoundaryPoints(polygon);
  const T interior = (area2 - boundary + 2) / 2;
  return interior;
}

////////////////////// point in polygon ////////////////////////////

template <class T>
int IntersectionType(const segment_t<T>& edge, const segment_t<T>& ray) {
  if (!Intersects(edge, ray)) return 0;
  return edge.p1.y > edge.p0.y ? +1 : -1;
}

template <class T>
int PointInsidePolygon(const polygon_t<T>& polygon, const point_t<T>& point) {
  vector<segment_t<T> > segments; int minX = polygon.P[0].x;
  for (int k = 0, N = polygon.N; k < N; ++k) {
    minX = min(minX, polygon.P[k].x);
    const segment_t<T> line(polygon.P[k], polygon.P[(k + 1) % N]);
    if (PointInSegment(line, point)) return 0;
    if (line.p0.y != line.p1.y) segments.push_back(line);
  }
  const segment_t<int> ray(point_t<int>(minX - 1, point.y), point);
  const int first = IntersectionType(segments[0], ray);
  int last = first, count = first != 0 ? 1 : 0;
  for (int k = 1, N = segments.size(); k < N; ++k) {
    const int next = IntersectionType(segments[k], ray);
    if (next != last && (last = next) != 0) ++count;
  }
  if (first == last && first != 0) --count;
  return count % 2 == 0 ? -1 : +1;
}

////////////////////// monotone chain convex hull ////////////////////////////

template <class T>
void EraseConcave(polygon_t<T>& hull, int fixed, const point_t<T>& point) {
  while (hull.N >= 2 && hull.N > fixed) {
    if (Orientation(hull.P[hull.N - 2], hull.P[hull.N - 1], point) > 0) break;
    --hull.N;
  }
}

template <class T>
int PushPoints(polygon_t<T>& hull, const polygon_t<T>& polygon, int fixed, int begin, int end) {
  const int dir = begin < end ? +1 : -1;
  for (int next = begin; next != end; next += dir) {
    EraseConcave(hull, fixed, polygon.P[next]);
    if (hull.N < 2 || Orientation(hull.P[hull.N - 2], hull.P[hull.N - 1], polygon.P[next]) > 0)
      hull.P[hull.N++] = polygon.P[next];
  }
  return hull.N;
}

template <class T>
polygon_t<T> GetConvexHull(polygon_t<T>& polygon) {
  sort(polygon.P, polygon.P + polygon.N, PointSort<T>());
  polygon.N = unique(polygon.P, polygon.P + polygon.N, PointEqual<T>()) - polygon.P;
  if (polygon.N == 1) return polygon; polygon_t<T> hull;
  PushPoints(hull, polygon, 0, 0, polygon.N);
  PushPoints(hull, polygon, hull.N, polygon.N - 2, 0);
  EraseConcave(hull, 0, polygon.P[0]);
  return hull;
}

////////////////////// voronoi diagrams ////////////////////////////

template <class T>
polygon_t<T> CutPolygon(const polygon_t<T>& polygon, const plane_t<T>& cutPlane) {
  polygon_t<T> cutPolygon; cutPolygon.N = 0;
  for (int k = 0; k < polygon.N; ++k) {
    const point_t<T>& p1 = polygon.P[(k + 0) % polygon.N];
    const point_t<T>& p2 = polygon.P[(k + 1) % polygon.N];
    const T dot1 = cutPlane.a * p1.x + cutPlane.b * p1.y + cutPlane.c;
    const T dot2 = cutPlane.a * p2.x + cutPlane.b * p2.y + cutPlane.c;
    if (dot1 * dot2 < T(0)) {
      const line_t<T> line = line_t<T>(p1, p2);
      const pair<int, T> cut = LinePlaneIntersection(cutPlane, line);
      const T t = cut.second;
      if (cut.first == POINT && t >= T(0) && t <= T(1)) {
        const point_t<T> cutPoint(p1.x + (p2.x - p1.x) * t, p1.y + (p2.y - p1.y) * t);
        cutPolygon.P[cutPolygon.N++] = cutPoint;
      }
    }
    if (dot2 >= T(0))
      cutPolygon.P[cutPolygon.N++] = p2;
  }
  return cutPolygon;
}

template <class T>
polygon_t<T> GetVoronoiPolygon(const vector<point_t<T> >& p, const polygon_t<T>& box, int k) {
  polygon_t<T> voronoi = box;
  for (int i = 0; i < p.size(); ++i) if (i != k) {
    const point_t<T> midPoint(T((p[k].x + p[i].x) / 2), T((p[k].y + p[i].y) / 2));
    const T a = T(p[k].x - p[i].x), b = T(p[k].y - p[i].y);
    const T c = -(a * midPoint.x + b * midPoint.y);
    const plane_t<T> cutPlane(a, b, c);
    voronoi = CutPolygon(voronoi, cutPlane);
  }
  return voronoi;
}

};

////////////////////// big integer basic operations ////////////////////////////

struct BigInteger {
  vector<int> digits;
  static const int BASE, DIGITS;
  BigInteger(int value = 0);
  BigInteger(const char* str);
  int get(int k) const;
  void set(int k, int value);
  void normalize();
  string toString() const;
};

const int BigInteger::BASE = 10000;
const int BigInteger::DIGITS = 4;

BigInteger::BigInteger(int value) {
  for (; value > 0; value /= BASE)
    digits.push_back(value % BASE);
}
BigInteger::BigInteger(const char* str) {
  char buffer[DIGITS + 1]; memset(buffer, '\0', sizeof(buffer));
  for (int length = strlen(str); length > 0; length -= DIGITS) {
    const int digitSize = min(DIGITS, length);
    memcpy(buffer, &str[length - digitSize], digitSize * sizeof(char));
    buffer[digitSize] = '\0';
    digits.push_back(atoi(buffer));
  }
  normalize();
}

int BigInteger::get(int k) const {
  if (k >= digits.size()) return 0;
  else return digits[k];
}
void BigInteger::set(int k, int value) {
  if (k >= digits.size()) digits.resize(k + 1, 0);
  digits[k] = value;
}
void BigInteger::normalize() {
  while (!digits.empty() && digits.back() == 0)
    digits.pop_back();
}

string BigInteger::toString() const {
  if (digits.empty()) return string("0", 1);
  char buffer[DIGITS + 1];
  const int length = sprintf(buffer, "%d", digits.back());
  string str(buffer, length);
  for (int k = digits.size() - 2; k >= 0; --k) {
    sprintf(buffer, "%0*d", DIGITS, digits[k]);
    str.append(buffer, DIGITS);
  }
  return str;
}

BigInteger operator <<(const BigInteger& b1, int offset) {
  BigInteger b2; b2.digits.insert(b2.digits.end(), offset, 0);
  b2.digits.insert(b2.digits.end(), b1.digits.begin(), b1.digits.end());
  b2.normalize(); return b2;
}

BigInteger& operator <<=(BigInteger& b1, int offset) {
  b1.digits.insert(b1.digits.begin(), offset, 0);
  b1.normalize(); return b1;
}

BigInteger operator >>(const BigInteger& b1, int offset) {
  BigInteger b2; if (offset < b1.digits.size())
    b2.digits.insert(b2.digits.begin(), b1.digits.begin() + offset, b1.digits.end());
  return b2;
}

BigInteger operator >>=(BigInteger& b1, int offset) {
  if (b1.digits.size() <= offset) b1.digits.clear();
  else b1.digits.erase(b1.digits.begin(), b1.digits.begin() + offset);
  return b1;
}

bool operator <(const BigInteger& b1, const BigInteger& b2) {
  if (b1.digits.size() < b2.digits.size()) return true;
  if (b2.digits.size() < b1.digits.size()) return false;
  for (int k = b1.digits.size() - 1; k >= 0; --k) {
    if (b1.digits[k] < b2.digits[k]) return true;
    if (b2.digits[k] < b1.digits[k]) return false;
  }
  return false;
}

bool operator >(const BigInteger& b1, const BigInteger& b2) {return b2 < b1;}
bool operator <=(const BigInteger& b1, const BigInteger& b2) {return !(b2 < b1);}
bool operator >=(const BigInteger& b1, const BigInteger& b2) {return !(b1 < b2);}
bool operator ==(const BigInteger& b1, const BigInteger& b2) {return b1.digits == b2.digits;}
bool operator !=(const BigInteger& b1, const BigInteger& b2) {return !(b1 == b2);}

BigInteger operator +(const BigInteger& b1, const BigInteger& b2) {
  BigInteger b3; int carry = 0;
  for (int k = 0; carry > 0 || k < b1.digits.size() || k < b2.digits.size(); ++k) {
    const int value = b1.get(k) + b2.get(k) + carry;
    b3.digits.push_back(value % BigInteger::BASE);
    carry = value / BigInteger::BASE;
  }
  return b3;
}

BigInteger& operator +=(BigInteger& b1, const BigInteger& b2) {
  int carry = 0;
  for (int k = 0; carry > 0 || k < b2.digits.size(); ++k) {
    const int value = b1.get(k) + b2.get(k) + carry;
    b1.set(k, value % BigInteger::BASE);
    carry = value / BigInteger::BASE;
  }
  return b1;
}

BigInteger operator -(const BigInteger& b1, const BigInteger& b2) {
  BigInteger b3; int borrow = 0;
  for (int k = 0; k < b1.digits.size() || k < b2.digits.size(); ++k) {
    const int value = BigInteger::BASE + b1.get(k) - b2.get(k) - borrow;
    b3.digits.push_back(value % BigInteger::BASE);
    borrow = 1 - value / BigInteger::BASE;
  }
  return b3.normalize(), b3;
}

BigInteger& operator -=(BigInteger& b1, const BigInteger& b2) {
  int borrow = 0;
  for (int k = 0; k < b1.digits.size() || k < b2.digits.size(); ++k) {
    const int value = BigInteger::BASE + b1.get(k) - b2.get(k) - borrow;
    b1.set(k, value % BigInteger::BASE);
    borrow = 1 - value / BigInteger::BASE;
  }
  return b1.normalize(), b1;
}

BigInteger operator *(const BigInteger& b1, int X) {
  BigInteger b2; int carry = 0;
  for (int k = 0; k < b1.digits.size(); ++k) {
    const int value = b1.digits[k] * X + carry;
    b2.digits.push_back(value % BigInteger::BASE);
    carry = value / BigInteger::BASE;
  }
  if (carry > 0) b2.digits.push_back(carry);
  return b2.normalize(), b2;
}

BigInteger& operator *=(BigInteger& b1, int X) {
  int carry = 0;
  for (int k = 0; k < b1.digits.size(); ++k) {
    const int value = b1.digits[k] * X + carry;
    b1.digits[k] = value % BigInteger::BASE;
    carry = value / BigInteger::BASE;
  }
  if (carry > 0) b1.digits.push_back(carry);
  return b1.normalize(), b1;
}

BigInteger operator *(const BigInteger& b1, const BigInteger& b2) {
  BigInteger b3;
  for (int k = 0; k < b2.digits.size(); ++k)
    b3 += (b1 * b2.digits[k]) << k;
  return b3;
}

BigInteger operator *=(BigInteger& b1, const BigInteger& b2) {
  BigInteger b3;
  for (int k = 0; k < b2.digits.size(); ++k)
    b3 += (b1 * b2.digits[k]) << k;
  return b1 = b3;
}

BigInteger operator /(const BigInteger& b1, int X) {
  assert(X > 0);
  BigInteger b2; b2.digits.resize(b1.digits.size());
  for (int k = b1.digits.size() - 1, carry = 0; k >= 0; --k) {
    carry = carry * BigInteger::BASE + b1.digits[k];
    b2.digits[k] = carry / X;
    carry %= X;
  }
  return b2.normalize(), b2;
}

BigInteger& operator /=(BigInteger& b1, int X) {
  assert(X > 0);
  for (int k = b1.digits.size() - 1, carry = 0; k >= 0; --k) {
    carry = carry * BigInteger::BASE + b1.digits[k];
    b1.digits[k] = carry / X;
    carry %= X;
  }
  return b1.normalize(), b1;
}

int operator %(const BigInteger& b1, int X) {
  assert(X > 0);
  int carry = 0;
  for (int k = b1.digits.size() - 1; k >= 0; --k)
    carry = (carry * BigInteger::BASE + b1.digits[k]) % X;
  return carry;
}

BigInteger& operator %=(BigInteger& b1, int X) {
  assert(X > 0);
  int carry = 0;
  for (int k = b1.digits.size() - 1; k >= 0; --k)
    carry = (carry * BigInteger::BASE + b1.digits[k]) % X;
  b1.digits.resize(1); b1.digits[0] = carry;
  return b1.normalize(), b1;
}

BigInteger operator /(const BigInteger& b1, const BigInteger& b2) {
  assert(b2.digits.size() > 0);
  BigInteger b3, carry; b3.digits.resize(b1.digits.size());
  for (int k = b1.digits.size() - 1; k >= 0; --k) {
    carry <<= 1; carry.set(0, b1.digits[k]);
    int lower = 0, upper = BigInteger::BASE;
    while (lower < upper) {
      const int middle = (lower + upper) / 2;
      if (b2 * middle <= carry) lower = middle + 1;
      else upper = middle;
    }
    b3.digits[k] = lower - 1;
    carry -= b2 * b3.digits[k];
  }
  return b3.normalize(), b3;
}

BigInteger operator /=(BigInteger& b1, const BigInteger& b2) {
  assert(b2.digits.size() > 0); BigInteger carry;
  for (int k = b1.digits.size() - 1; k >= 0; --k) {
    carry <<= 1; carry.set(0, b1.digits[k]);
    int lower = 0, upper = BigInteger::BASE;
    while (lower < upper) {
      const int middle = (lower + upper) / 2;
      if (b2 * middle <= carry) lower = middle + 1;
      else upper = middle;
    }
    b1.digits[k] = lower - 1;
    carry -= b2 * b1.digits[k];
  }
  return b1.normalize(), b1;
}

BigInteger operator %(const BigInteger& b1, const BigInteger& b2) {
  assert(b2.digits.size() > 0); BigInteger carry;
  for (int k = b1.digits.size() - 1; k >= 0; --k) {
    carry <<= 1; carry.set(0, b1.digits[k]);
    int lower = 0, upper = BigInteger::BASE;
    while (lower < upper) {
      const int middle = (lower + upper) / 2;
      if (b2 * middle <= carry) lower = middle + 1;
      else upper = middle;
    }
    carry -= b2 * (lower - 1);
  }
  return carry;
}

BigInteger operator %=(BigInteger& b1, const BigInteger& b2) {
  assert(b2.digits.size() > 0);
  BigInteger carry;
  for (int k = b1.digits.size() - 1; k >= 0; --k) {
    carry <<= 1; carry.set(0, b1.digits[k]);
    int lower = 0, upper = BigInteger::BASE;
    while (lower < upper) {
      const int middle = (lower + upper) / 2;
      if (b2 * middle <= carry) lower = middle + 1;
      else upper = middle;
    }
    carry -= b2 * (lower - 1);
  }
  return b1 = carry, b1;
}

//////////////////////////////////// big integer advanced functions //////////////////////////////////////////////////////

BigInteger Power(const BigInteger& b, int p) {
  if (p == 0) return BigInteger(1);
  if (p == 1) return b;
  BigInteger x = Power(b, p / 2); x *= x;
  if (p % 2 != 0) x *= b; return x;
}

int GetNthRootDigit(const BigInteger& X, BigInteger& Y, int N) {
  int lower = 0, upper = BigInteger::BASE;
  while (lower < upper) {
    const int middle = (lower + upper) / 2;
    Y.set(0, middle);
    if (Power(Y, N) <= X) lower = middle + 1;
    else upper = middle;
  }
  return lower - 1;
}

BigInteger NthRoot(const BigInteger& b1, int N) {
  BigInteger X(0), Y(0), R(0);
  for (int B = (b1.digits.size() + N - 1) / N - 1; B >= 0; --B) {
    X <<= N, Y <<= 1;
    for (int k = 0; k < N; ++k)
      X.set(k, b1.get(B * N + k));
    X.normalize();
    Y.set(0, GetNthRootDigit(X, Y, N));
    R = X - Power(Y, N);
  }
  return Y;
}

//////////////////////////////////// integer fractions //////////////////////////////////////////////////////

template <class T> struct fraction_t {
  T N, D;
  fraction_t() : N(0), D(1) {}
  fraction_t(T N_) : N(N_), D(1) {}
  fraction_t(T N_, T D_) : N(N_), D(D_) {normalize();}
  void normalize() {DivByGCD(N, D); if (D < 0) N = -N, D = -D;}
};

template <class T> fraction_t<T> operator +(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G = GCD(f1.D, f2.D), D = (f1.D / G) * f2.D;
  const T N = f1.N * (f2.D / G) + f2.N * (f1.D / G);
  return fraction_t<T>(N, D);
}
template <class T> fraction_t<T> operator -(const fraction_t<T>& f) {
  return fraction_t<T>(-f.N, f.D);
}
template <class T> fraction_t<T> operator -(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G = GCD(f1.D, f2.D), D = (f1.D / G) * f2.D;
  const T N = f1.N * (f2.D / G) - f2.N * (f1.D / G);
  return fraction_t<T>(N, D);
}
template <class T> fraction_t<T> operator *(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G1 = GCD(f1.N, f2.D), G2 = GCD(f2.N, f1.D);
  const T N = (f1.N / G1) * (f2.N / G2);
  const T D = (f2.D / G1) * (f1.D / G2);
  return fraction_t<T>(N, D);
}
template <class T> fraction_t<T> operator /(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  assert(f2.N != 0);
  const T G1 = GCD(f1.N, f2.N), G2 = GCD(f1.D, f2.D);
  const T N = (f1.N / G1) * (f2.D / G2);
  const T D = (f2.N / G1) * (f1.D / G2);
  return fraction_t<T>(N, D);
}
template <class T> fraction_t<T>& operator +=(fraction_t<T>& f1, const fraction_t<T>& f2) {return f1 = f1 + f2;}
template <class T> fraction_t<T>& operator -=(fraction_t<T>& f1, const fraction_t<T>& f2) {return f1 = f1 - f2;}
template <class T> fraction_t<T>& operator *=(fraction_t<T>& f1, const fraction_t<T>& f2) {return f1 = f1 * f2;}
template <class T> fraction_t<T>& operator /=(fraction_t<T>& f1, const fraction_t<T>& f2) {return f1 = f1 / f2;}

template <class T> bool operator <(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G1 = GCD(f1.N, f2.N), G2 = GCD(f1.D, f2.D);
  return (f1.N / G1) * (f2.D / G2) < (f2.N / G1) * (f1.D / G2);
}
template <class T> bool operator >(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G1 = GCD(f1.N, f2.N), G2 = GCD(f1.D, f2.D);
  return (f1.N / G1) * (f2.D / G2) > (f2.N / G1) * (f1.D / G2);
}
template <class T> bool operator <=(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G1 = GCD(f1.N, f2.N), G2 = GCD(f1.D, f2.D);
  return (f1.N / G1) * (f2.D / G2) <= (f2.N / G1) * (f1.D / G2);
}
template <class T> bool operator >=(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G1 = GCD(f1.N, f2.N), G2 = GCD(f1.D, f2.D);
  return (f1.N / G1) * (f2.D / G2) >= (f2.N / G1) * (f1.D / G2);
}
template <class T> bool operator ==(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G1 = GCD(f1.N, f2.N), G2 = GCD(f1.D, f2.D);
  return (f1.N / G1) * (f2.D / G2) == (f2.N / G1) * (f1.D / G2);
}
template <class T> bool operator !=(const fraction_t<T>& f1, const fraction_t<T>& f2) {
  const T G1 = GCD(f1.N, f2.N), G2 = GCD(f1.D, f2.D);
  return (f1.N / G1) * (f2.D / G2) != (f2.N / G1) * (f1.D / G2);
}

};

int main() {
  return 0;
}
